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Discussion 


The  p-version  of  the  finite  element  method  is  a  new  approach  to  finite 
element  analysis  which  has  been  demonstrated  to  lead  to  significant  computa¬ 
tional  savings,  often  by  orders  of  magnitude  (This  approach  was  formerly 
called  the  constraint  method;  the  new  term  p-version  is  more  descriptive). 
Conventional  approaches  (called  the  h-version)  generally  employ  low  order 
polynomials  as  basis  functions.  Accuracy  is  achieved  by  suitably  refining 
the  approximating  mesh.  The  p-version  uses  polynomials  of  arbitrary  order 
p  >  2  for  problems  in  plane  elasticity  where  CO  continuity  is  required  and 
polynomials  of  order  p  >  5  for  problems  in  plate  bending  where  Cl  continuity 
is  required. 

Hierarchic  elements  which  implement  the  p-version  efficiently  are  used 
together  with  precomputed  arrays  of  elemental  stiffness  matrices. 

Major  accomplishments  during  this  past  year  are  summarized  in  the  follow¬ 
ing  three  documents  which  are  enclosed: 

1.  "Comparative  Rates  of  h-  and  p-  convergence  in  the  Finite  element  Anal¬ 
ysis  of  a  Model  Bar  Problem"  by  I.  Norman  Katz  (abstract),  presented 

at  SIAM  1978  Fall  Meeting,  October  30,  31,  November  1,  1978  in  Knoxville, 
Tennessee. 

2.  "The  p-Version  of  the  Finite  Element  Method"  by  I.  Babuska,  B.  S.  Szabo, 
and  I.  Norman  Katz  (Report,  submitted  for  publication).  Report  WU/CCM-79/1, 
May  1979 

3.  "Hierarchic  Families  for  the  p-Version  of  the  Finite  Element  Method",  by 
I.  Babuska,  I.  Norman  Katz  and  B.  A.  Szabo,  Proceedings  of  the  Third 
IMACS  International  Symposium  on  Computer  Methods  for  Partial  Differen¬ 
tial  Equations,  Lehigh  University,  Bethlehem,  PA,  June  20  -  22,  1979 
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I .  Norman  Katz ,  Washington  University,  St.  Louis, 
Missouri  63130 


Comparative  Rates  of  h-  and  p-convergence  in  the 
Finite  Element  Analysis  of  a  Model  Bar 
Problem 

The  conventional  approach  to  finite  element  stress 
analysis  of  a  body  defined  by  a  polygonal  domain 
ft  (in  two  dimensions)  is  to  triangulate  and  to 
seek  accuracy  by  letting  h,  the  maximum  diameter 
of  all  elements  in  the  triangulation,  tend  to 
zero.  This  approach,  called  h-convergence ,  has 
been  the  subject  of  intensive  investigation. 
Another  approach  which  is  being  developed  at 
Washington  University  is  to  .fix  the  triangulation 
of  (2  and  to  let  p,  the  degree  of  the  complete, 
conforming,  approximating  polynomial  over  each 
triangle,  tend  to  infinity.  Extensive  numerical 
tests  have  shown  that  the  second  approach,  called 
p- convergence ,  is  considerably  more  accurate  than 
the  first,  even  in  problems  whose  solutions  have 
singularities  such  as  cracks  or  corners. 

In  order  to  illustrate  the  comparative  rates  of 
convergence,  a  model  (one-dimensional)  bar  prob¬ 
lem  is  studied.  Asymptotic  analysis  leads  to  ex¬ 
pressions  for  the  rates  of  convergence  in  the  two 
approaches,  when  the  solution  possesses  a  singu¬ 
larity  which  is  known  a  priori.  It  is  demon¬ 
strated  that  the  order  of  p-convergence  is  twice 
that  of  h-convergence,  provided  that  the  singular¬ 
ity  is  located  at  some  node  of  a  finite  element. 
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ABSTRACT 

In  the  p-version  of  the  finite  element  method  the  triangulation  is 
fixed  and  the  degree  p  of  the  piecewise  polynomial  approximation  is 
progressively  increased  until  some  desired  level  of  precision  is  reached. 

In  this  paper  we  first  establish  the  basic  approximation  properties 
of  some  spaces  of  piecewise  polynomials  defined  on  a  finite  element 
triangulation.  These  properties  lead  to  an  a  priori  estimate  of  the 
asymptotic  rate  of  convergence  of  the  p-version.  The  estimate  shows  that 
the  p-version  gives  results  which  are  not  worse  than  those  obtained  by  the 
conventional  finite  element  method  (called  the  h-version,  in  which  h  rep¬ 
resents  the  maximum  diameter  of  the  elements)  when  quasi-uniform  triangula¬ 
tions  are  employed  and  the  basis  for  comparison  is  the  number  of  degrees 
of  freedom.  Furthermore,  in  the  case  of  a  singularity  problem  we  show 
(under  conditions  which  are  usually  satisfied  in  practice)  that  the  rate  of 
convergence  of  the  p-version  is  twice  that  of  the  h-version  with  quasi-uniform 
mesh.  Inverse  approximation  theorems  which  determine  the  smoothness  of  a 
function  based  on  the  rate  at  which  it  is  approximated  by  piecewise  poly¬ 
nomials  over  a  fixed  triangulation  are  proved  both  for  singular  and  non¬ 
singular  problems. 

We  present  numerical  examples  which  illustrate  the  effectiveness  of 
the  p-version  for  a  simple  one  dimensional  problem  and  for  two  problems  in 
two  dimensional  elasticity.  We  also  discuss  round  off  error  and  computa¬ 
tional  costs  associated  with  the  p-version.  Finally  we  describe  some 
important  features,  such  as  hierarchic  basis  functions,  which  have  been 
utilized  in  COMET-X,  an  experimental  computer  implementation  of  the 
p-version. 
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1. 


INTRODUCTION 


The  finite  element  method,  one  of  the  most  widely  used  numerical 
methods  for  solving  certain  types  of  differential  equations,  is  based 
on  approximating  the  solution  by  piecewise  smooth  functions,  speci¬ 
fically  polynomials,  on  convex  subdomains  such  as  triangles.  In  general, 
the  degree  of  the  polynomials  is  fixed  at  some  arbitrarily  chosen  low 
number.  No  consensus  exists  at  the  present  time  concerning  the  most 
suitable  (optimal)  degree  p  of  the  polynomials. 

The  mathematical  justification  of  the  finite  element  method  is 
based  on  asymptotic  analyses  in  which  p  is  kept  bounded  and  the  diameters 
of  the  element  subdomains  approach  zero.  However,  it  has  been  observed 
by  several  investigators  that  the  sizes  of  elements  used  in  practical 
computations  are  often  outside  of  the  range  of  asymptotic  behavior. 

Because  the  maximum  diameter  of  finite  elements  is  usually  denoted 
by  h,  we  shall  refer  to  this  (conventional)  approach  as  the  h- vers ion 
of  the  finite  element  method. 

From  the  theoretical  point  of  view  one  can  justify  the  finite 
element  method,  also  in  the  asymptotic  sense,  when  the  subdomains  are 
kept  constant  and  the  degree  of  the  approximating  polynomials  tends  to 
infinity.  We  shall  refer  to  this  method  of  approximation  as  the  p-version 
of  the  finite  element  method. 

The  p-version  of  the  finite  element  method  is  similar  to  the  Ritz 
method  but  there  is  one  very  important  difference:  In  the  p-version  of 
the  finite  element  method  the  domain  of  interest  is  divided  into  convex 
subdomains  and  the  polynomial  approximants  are  piecewise  smooth  only 
over  individual  convex  subdomains.  In  the  Ritz  method,  on  the  other 
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hand,  the  solution  over  the  entire  domain  is  approximated  by  smooth 
functions.  This  difference  accounts  for  the  greater  versatility  and 
higher  rate  of  convergence  of  the  p-version  of  the  finite  element 
method  over  both  the  Ritz  method,  and  the  h-version  of  the  finite 
element  method,  as  demonstrated  here. 

In  this  paper  we  analyze  the  p-version  of  the  finite  element 
method  and  its  theory,  and  discuss  the  implementation  characteristics 
of  the  method  based  on  the  computer  program  COMET-X,  developed  during 
the  last  few  years  at  Washington  University  in  St.  Louis.  We  also 
examine  the  potential  for  further  development  of  the  p-version.  We 
remark  that,  from  the  computational  point  of  view,  and  from  the  point  of 
view  of  the  architecture  of  the  computer  program,  there  are  significant 
differences  between  the  p-version  when  p  is  in  the  range  of  6,7,8  and  the 
h-version  when  p  is  in  the  range  1,2,3. 

We  present  a  proof  for  the  rate  of  convergence  in  the  p-version 
and  show  that  the  polynomials  are  able  to  "absorb"  singularities, 

including  e.g.,  corner  singularities,  when  they  are  located  at  the  vertices 
of  triangles.  This  does  not  occur  when  the  corner  singularities  are 
not  located  at  vertices. 

Comparison  of  the  asymptotic  behavior  of  the  h-version,  based  on 
uniform  or  quasi-uniform  mesh  refinement  on  one  hand,  and  the  p-version 
on  the  other,  the  basis  of  comparison  being  the  number  of  degrees  of 
freedom,  shows  that  the  r;  of  convergence  of  the  p-version  cannot  be 
slower  than  the  rate  of  convergence  of  the  h-version  and,  furthermore, 
when  corner  singularities  are  present  at  vertices,  the  rate  of  convergence 
of  the  p-version  is  exactly  twice  that  of  the  h-version. 
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2.  BASIC  NOTATIONS 

2 

Throughout  this  paper  R  will  be  the  two  dimensional  Euclidean 
2  2 

space  (x^,X2)  =  x  E  R  ,  OCR”  will  be  a  bounded  domain  with  piecewise 
smooth  boundary  30.  In  particular  we  will  deal  with  polygonal  domains 
(We  exclude  -  for  technical  reasons  -  the  slit  domain,  although  the 
results  of  this  paper  can  be  generalized  to  this  case  too  with  some,  but 
not  essential,  technical  difficulties). 

E(0)  shall  be  the  space  of  all  real  C  functions  on  0,  with  continuous 
extensions  of  all  derivatives  on  0.  'All  functions  of  E(0)  with  compact 
support  in  0  form  a  subspace  P(I2)  C  E(fi).  As  usual,  L2  (f2)  =  will 

be  the  space  of  all  square-lntegrable  functions  on  Q  with  the  inner 
product . 


(u,v) 


o,n 


uvdx  y 


dx  =  dx^dx^ 


and  the  corresponding  norm  ||*||  .  In  addition  for  any  k  >  1,  integral, 

Ir  Ir  _ 

the  Sobolev  spaces  H  (fl)  resp  Hq(Q)  will  be  the  completions  of  E(K) 
resp.  under  the  norm 


where 


'k.n 


E 

0< la  I <k 


“l  “2 
3x^  8X2 


a  =  (a1,a2) 
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>  0  integral  i  *  1,2  and  |  ot  |  =  +  a^.  The  inner  product  in  H^(SJ) 

J, 

will  be  denoted  by  .  For  k  >  0  nonintegral  the  spaces  H  (Q) 

and  Hg(ft)  are  defined  by  usual  interpolation.  More  precisely  for  k  =  k  +  9, 

k 

0  <  9  <  1,  we  define  =  [H  h'<0+'*']q  .  by  application  of  the  usual 

o ,  z 

K-method  of  interpolation  (For  more  see  [7]).  [Other  notations  are 
k  k  k 

H  =  2  w^ere  2  t*ie  usua^  Besov  space] . 

For  p  >  0  we  write 

Q(p)  =  {x1,x2|  | x1 |  <  p,  | X^ |  <  P>  , 

Q(p)  =  {x1,x2|0  <  x1  <  p,  0  <  x2  <  p} 

and  by  EpER(Q(p))C  E(Q(p))  we  denote  the  space  of  all  functions  with 
period  2p  and  by  H^RR(Q(p))  its  closure  in  Hk(Q(p)). 

We  will  deal  also  with  Sobolev  spaces  in  one  dimension.  Analogously 
as  before  we  will  denote  by 


I(P) 


1*3.1  <  P> 


and  Hk(I(p)),  Hg(I(p)),  H^ER(I(p))  will  have  the  obvious  meaning. 

Finally  we  need  to  introduce  the  spaces  C  E(ft)  of  all  algebraic 

polynomials  of  degree  not  higher  than  p  and  F^(Q(p))  (resp.  F^(I(p))  the 
space  of  all  trigonometric  polynomials  of  degree  at  most  p  and  period  2p. 


3. 


THE  CONCEPT  OF  P-CONVERGENCE  OF  THE  FINITE  ELEMENT  METHOD 


3.1  The  model  problem 

We  will  be  interested  in  the  model  problem 

-Au  +  U  *=  f  on  n0,  f  e  H°(n0),  (3.1) 


Fu  *  o  on  3nQ  ,  (3.2) 

where  is  a  bounded  polygonal  domain  and  Fu  =  u  or  Tu  *  We  can 

u  3n 

easily  generalize  our  results  also  to  other  boundary  conditions.  As 
usual  we  will  interpret  the  problem  (3.1),  (3.2)  in  a  weak  sense,  namely 
we  seek  uQ  E  hq(^q)  resp.  U()  6  H1  (^Q)  so  that 


B(Vv)  =  (f’V)0,D0 

for  all  ve  hJ(S20)  resp.  v  E  H1(nQ) 
where  we  have  used  the  notation 


(3.3) 


B(VV>  .  <»0,»)liS, 


0  - 


(3.4) 


Uq  satisfying  (3.3)  obviously  exists  and  is  uniquely  determined. 


3.2  Description  of  the  p-version  of  the  finite  element  method 

Let  S  be  a  (fixed)  triangularization  of  JJ  ,  S  -  {T  },  i  -  l,...,m 


where  T^  are  (open)  triangles  such  that  U  T±  «  and  T_, 


i-  1  *  i 
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have  either  a  common  (entire)  side  or  a  vertex  or  T  Pi  T.  =  d>.  Denote 

i  J 

now  by  Pp^  (Qq)  C  H^(Qg)  the  subset  of  all  functions  u€  H^(Dq)  such 

that  if  u N  is  the  restriction  of  u  on  T.,  then  we  have  u,_  .  £  P  (T.) 

(T.)  i’  (T  )  pv  i' 

[S]  1  1 

i.e.  Pp  J (Dq)  consists  of  all  functions  which  are  piecewise  polynomials 

1  fSI 

of  degree  at  most  p  and  belong  to  H  (OJ .  Further  let  Pl  * 

u  p,  u  o 

ppS1<Vn  HJ‘V- 


The  concept  of  the  p-version  of  the  finite  element  method  consists  of 

finding  u  p  =  1,2,....  u  E  P^^(£lrt)  (resp.  P^  (ft,,))  (for  the  boundary 
p  pp,UU  pU 

condition  Tu  =  0  resp.  Tu  =  -|^)  so  that  (3.3)  holds  for  any  v  G  P^|(J)„) 

dn  p,U  U 

(resp.  PpS^(QQ)). 

Study  of  the  p-version  of  the  finite  element  method  was  initiated 
at  the  School  of  Engineering  and  Applied  Science  of  Washington  University 
in  St.  Louis  [25] in  1970.  It  has  been  implemented  there  in  various 
aspects  of  stress  analysis  with  very  good  results,  particularly  in 
connection  with  linear  elastic  fracture  mechanics.  Development  of  the 
p-version  is  continuing  at  the  Center  for  Computational  Mechanics  at 
Washington  University. 


3.3  The  basic  approximation  properties  of  Pp^  (D^)  and  Pp^g(%)  ‘ 

THEOREM  3.1.  Let  u  E  Hk(f2-).  Then  there  exists  a  sequence  z  E  P^^  (ft.)  , 

U  p  p  0 

p  =  1,2...  such  that  for  any  0  <  £  <  k  (£,k  not  necessarily  integral) 


u-z 


P 


<  cP‘(k'-e) 


(3.5) 


where  C  is  independent  of  u  and  p  (it  depends  e.g.  on  £  and  k  etc.). 

Proof.  The  proof  is  a  standard  one.  First  we  prove  it  for  £  and  k 
integral.  We  will  construct  zp  E  Pp(nQ)  such  that  (3.5)  is  satisfied. 


I 


I 


Assume  that  QQ  C  Q(pq).  Because  is  a  polygon,  it  is  a  Lipschitzian 
domain  and  therefore  there  exists  an  extension  U  of  u  E  Hk(nQ)  on  Q(2pQ) 
such  that  supp  U  C  Q ( 3/ 2  pQ)  and 


|u|  lk,Q(2p0)  -  CI  'u'  'k.r.Q 


with  C  independent  of  u.  As  usual  we  have  U  =  Xu  where  T  is  a  linear 
mapping  of  H0(flQ)  into  H°(Q(2pQ)),  (see  e.g.  [24])  (which  also  maps 
H°(Qq)  into  H°(Q(2p0)). 

Now  let  $  be  the  (one  to  one)  mapping  of  Q(tt/2)  onto  Q(2p^  )  determined 
by  the  transformation  of  coordinates  £  =  (^.S^  £  Q(*/2),  x  =  (x^.x^  E  Q(2pQ) 


xi  =  2p0  Sin  5i  1  =  1,2 


written  in  the  form 


$(?)  =  x. 


V(5)  -  0(*«)) 


and  let 


Q  »  *l"iJ[Q(3/2  pQ)]  C  Q(ir/2) 


(il  is  the  inverse  mapping  to  $) ,  We  have 


I 
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Supp  VC  Q. 


Obviously  the  mapping  $  is  a  regular  analytic  one-to-one  mapping  of  Q 
onto  Q(3/2  Pq).  Now  V  £  Hp^CQCir)),  is  symmetric  with  respect  to  the 
lines  5  *  ±  ir/2  and  using  (3.6)  and  (3.9)  we  obtain 


>k,Q(ir)  -  C  *  'k,C0  ' 


(3.10) 


It  is  well  known  that  the  partial  sum  t  of  the  Fourier  series  of  V 

P 

gives  the  sequence  of  trigonometric  polynomials  t  £  P  (Q(ir))  such  that 

P  P 


for  k  >  t 


iv-'p|  i'p'<k'wnviik><jW 


icp-<k-«iNI  k_a. 


(3.11) 


tp  are  obviously  symmetric  with  respect  to  the  lines  =  ±  tt/2  as  V  is. 

It  is  readily  seen  that  tp(£)  =  zp($(£)),  where  zp  is  an  algebraic  polynomial 

of  degree  not  higher  than  p.  Because  *  is  a  regular,  analytic  mapping 

of  Q  onto  Q (3/ 2  PQ)  (3.11)  yields  (3.5)  for  k,£  integral. 

Now  let  us  generalize  our  result  to  t,k  not  integral.  Recall  that 

for  given  (fixed)  p  the  polynomial  z  was  constructed  from  a  linear  map 

P 

L  ,  L  u  “  z  ,  where  L  is  a  linear  mapping  of  H^(n„)  into  P  (ft„)  satisfying 

p  p  p  p  0  p  0 

(3.5)  for  £,k  integral.  Applying  general  interpolation  theory  we  get 
(3.5)  for  all  0  <  l  <  k. 


The  proof  of  the  next  theorem  is  more  complicated. 
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THEOREM  3.2.  Let  u£  H^(3q)  D  Hq(J2q).  There  exists  a  sequence  z^£  P^q(Qq) 
P  =  1,2,...  such  that  for  any  k  >  1  (not  necessarily  integral)  and  any 
£  >  0 


ll  $1  <  °P 


-(k-l)+€i 


(3.12) 


where  C  is  independent  of  p  and  u  (it  depends  e.g.  on  eand  k) . 

Remark  1.  In  contrast  to  Theorem  3.1  the  statement  is  false  if  only 

r  q  ■) 

P  0(^0)  instead  of  PL  i(nn)  is  considered.  This  is  easy  to  see  if 
p,u  U  p,U  0 

Hq  is  e.g.  an  L-shaped  domain  as  shown  in  Fig.  3.1. 


Figure  3.1 
An  L-shaped  domain 

In  fact  any  u  E  Pp  0(ftQ)  is  zero  in  (x^.O),  0  <  x1  <  1  and  therefore  - 
because  it  is  a  polynomial  -  has  to  be  zero  on  the  entire  line  (x^O), 

-1  <  x^  <  1.  This  of  course  leads  to  a  contradiction  because  of  Sobolev's 
imbedding  theorem  of  Hg(flQ)  into  if  (1(1))  . 
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Remark  2.  It  is  not  clear  whether  the  terms  in  (3.12)  can  be  removed. 
Remark  3.  The  theorem  can  be  stated  more  generally.  We  have  restricted 
ourselves  to  this  case  (i.e.  ]|  •  ||^  ^)  only  because  it  is  sufficient  for 
our  purpose. 

Before  proving  theorem  3.2  we  will  state  a  lemma. 

Lemma  3.1.  Let  S  be  a  triangle  with  vertices  ,  and  sides  s.,  i  =  1,2,3 
(see  Fig.  3.2) 


A  3 


Figure  3.2 

The  typical  general  triangle 

Let  v  E  P  (s  ).  Then  there  exists  V  £  P  (S)  such  that  V  =  0  on  s _  and 

p ,  u  x  p  2 

V  =  v  on  and 

I  M  I i,s  -  c  Mvi  ll,  s  (3.13) 

where  C  (dependent  on  S)  is  independent  of  v  and  p. 

Remark.  By  v  G  P^  g(s^)  we  mean  of  course  a  polynomial  in  the  variable 
s2  so  that  v  *  0  at  the  end  points  of  s^,  the  vertices  A2,  A^. 
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Proof.  Without  any  loss  of  generality  we  can  assume  that  S  is  the 
triangle  shown  in  Fig.  3.3  with  vertices  (0,0)  (1,0)  (1,1) 


Figure  3.3 

The  Standard  Triangle 

Then  s^  =  (x^,0),  0  <  x^  <  1.  Because  v(x)  is  a  polynomial  and  because 
v(0)  =  v(l)  =  0  by  assumption  we  have 


v(x1)  =  x1(l-x1)v1(x1) . 


with  v1(x1)  a  polynoimal  of  degree  at  most  p-2.  Define 

(x,-x  ) 

V(x)  -  V(x1(x2)  -  v(xx)  — i — — 


(3.14) 


Obviously  V  £  P  (S)  ,  V  =  0  on  s2J  s^, 
X2 

—  is  bounded  on  S  we  get  | |v| |Q  s  < 


and  V  =  v 


on 


s 

1’ 


Finally  because 


Writing 
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^<f-)2  -fo  dx2  f  (f-)2  dx1 


2  2 

we  easily  get  (3.13)  when  using  the  obvious  inequality  v  (x  )  <  x  | |v| |  s 

1  1  ^  *  1 

and  the  lemma  is  completely  proven. 

Now  we  can  prove  Theorem  3.2. 

Proof  of  Theorem  3.2.  1)  Let  k  >  2.  Applying  theorem  3.1  there 

exists  z  E  P  (Q„)  such  that 
P  P  0' 


lu"2plU>flo  -  Cp_(k"£)l  lul  lk,n0’  k 


(3.15) 


r  i  ] 

Let  {xLJ  }  j  =  l,...m  be  the  set  of  all  vertices  of  the  triangles  T^E  S 

1+E, 

belonging  to  3S1.  Because  H  (.Jig)  >  e  >  0  is  imbedded  in  the  space  of 
continuous  functions,  we  can  obviously  modify  z^  to  z*  by  subtracting  a 
polynomial  of  fixed  degree  p^  <  m  so  that 


I  *ii  -(k-£),  — (k— 1— e).  M  , 

lu_zpl  -  c(p  p  '  >llul 


(3.16) 


with  e  >  0  arbitrary.  In  fact  z  =  z^-z*  a  polynomial  of  fixed  (inde¬ 
pendent  of  p)  degree  m  determined  by  its  values  at  the  points  (x^  }. 

By  theorem  3.1  we  have  |z  (x  ^  ^)  |  <  Cp~^_1^+£|  |u|  |  and  so 

P  ~ 

i  I**,  i  -  He- 1)4*  e  a  u 

\\z\\p  o  <  Cp  for  all  0  <  &  <  m.  Because  obviously 

„  ’  0  ~ 

I  I z 1 1  o  =  |  | z  j  J  0  for  any  r  >  m,  (3.16)  follows  readily. 

r,>4g  m.iig 

We  see  that  on  every  side  s  C  3Ci  of  T  E  S  we  have  u  =  0  and  there¬ 
fore  ! I zp llis  -  | I u— zp I  I  2  fl  applying  the  Sobolev  imbedding  theorem. 

Using  lemma  3.1  we  can  now  find  z**  E  P^  (fl_)  so  that  z°  *  z*  -  z**  E  P^I(ft_) 

P  P  0'  P  p  P  rp,0v“0/ 

and 
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u-z 


lia  <c| 
o 


—  2* 


U-2 


l2  a  <  °p 


-a-2) 


'k,nr 


(3.17) 


where  C  depends  on  k  >  2  but  is  independent  of  u.  (3.17)  can  obviously 
be  written  as 


u-z 


p1  li,n 


o 


<  Cp 


.(k-i)(i-^) 


(3.18) 


2)  Let  now  be  the  orthogonal  projection  in  the  scalar  product  of  H1 
of 


“k<v n  »£<v  p$<%> 


Let  z  =  R  u. 
P  P 


We  obviously  have 


z  -u 
P 


i,n 


o 


< 


(3.19) 


and  from  (3.18)  for  k  >  2  we  have 


<  C (k)p 


-(k-1)  (1  -  ^) 


(3.20) 


For  1  <  s  <  k,  let 


HS(G0) 


[hJ,  Hk«i0)  n  hJ] 


S  — 1 

k-1 


2 


We  obtain  by  applying  interpolation  theory 
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|z_-u| |,  „  <  C(k,s)  p  Uj | u | 


'P  'i,nr 


HS(nQ) 


with 


(3.22) 


u  =  (k-1)  (1  -  j^j)(~)  =  (s-l)Cl  -  j^)  .  (3.23) 

Therefore  given  e  >  0  and  s  >  1  we  can  select  k^  so  that 

(s-i)(i  -k4r)  >  (s-D  -  e  *  0.24) 

and  so  (3.12)  holds  when  the  norm  ||u||_,  ,  instead  of  ||u||  ,  is 

Hk(nQ)  k,n0 

used. 

On  the  other  hand  from  [3],  see  also  [29],  it  follows  that  when  Qq 

~k  k  1 

is  a  polygon, then  the  spaces  H  (0^)  and  H  (fi^)  D  are  equivalent. 

This  completes  the  proof. 

3.4  The  inverse  approximation  theorem 

We  have  proven  theorems  about  approximability  properties  of  the 
spaces 


■fhv  “<*  ppfi<V  • 


Now  we  will  prove  the  inverse  approximation  theorem. 

THEOREM  3.3.  Let  vE  H  (Q(p))  and  let  there  exists  a  sequence  of  poly¬ 
nomials  z  E  P  (Q(p)),  p  *  1,2,...  such  that 
P  P 

lk.QCp)  i  ^7 


r  >  0 


(3.25) 
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k  >  0  integral.  Then  v  S  H  (Q(p*)),  p*  <  p,  e  >  0  arbitrary  and 
after  restriction  of  v  onto  Q(p*) 


|Vi ik+r-€,Q(P*)  -  A(p,p*,k,r,E)  fiivil0>Q(p)  +  K1  (3.26) 


Proof.  Let  cu  =  (x  -p“)  (x?-p-)  Then  writing  v*  =  vco,  z*  =  z  oj 

L  p+4(k+2)  p 


V*(5)  =  v*(KO) 


'p+4  (k+2)  ^  Zp+4(k+2)  e  Fp+4(k+2) 


$(C)  =  x 


x.  =  p  sxn 

l  i 


we  obtain 


p+4k+2vs/  I  lk>Q(ir)  1  ^ 


(3.27) 


Now  using  theorem  5.4.1  p.  200  of  [17]  [For  a  proof  using  the  basic 
interpolation  theory  directly  see  e.g.  [3]]  it  follows  that 


^ k+r-c ,Q(ir)  -  A[l  lv*l  Ifc.Q^)  +  K] 
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(More  precisely  by  the  mentioned  theorems  we  obtain  the  norm  of  V*  in 
the  Nikolsky  spaces  B  +r(Q(ir)  which  majorizes  the  norm  ||  ||  .  .  .  )  . 

Now  using  the  fact  that  the  mapping  $  is  a  one-to-one  analytic 
one  on  Q(P*)  and  w(x)  >  a  >  0  on  Q(p*)  we  immediately  obtain  (3.26). 
Inequality  (3.25)  holds  only  on  Q(p*)  and  in  general  is  not  true  on  Q(p). 
Nevertheless  we  can  prove  the  next  theorem. 

THEOREM  3.4.  Let  v  £  H  (Q(p))  and  suppose  that  the  other  assumptions  of 
Theorem  3.3  are  satisfied,  then  v£  yk+r/2  e(Q(p))  an(j 

^V^k+r/2-£,Q(p)  -  A(E)[I  lvl  lkjQ(p)  +  K]  (3.28) 

The  proof  of  this  theorem  is  a  consequence  of  the  above  mentioned 
theorem  5.4.1  in  [17]  provided  that  for  integral  k  >  0  the  following 
inequality  of  Bernstein  type 


P1 'k.Q(p) 


„  2k 

<  Cp  M  z 


P1  0,Q(p) 


(3.29) 


holds  for  any  z  £  P  (Q(p))  with  C  independent  of  p  and  z  . 

P  P  p 

Let  us  remark  that  in  the  case  of  trigonometrical  polynomials  we 

k  2k 

have  in  (3.29)  the  term  p  instead  p  .  We  will  prove  (3.29)  in  the 
next  few  lemmas. 

Lemma  3.2.  Let  z  (x) ,  x  £  1(1)  be  a  polynomial  (in  one  variable)  of 
P 

degree  p.  Then 


1 

,s 
d  z 

1 

P 

dx 

s 

0,i  1  C(s)p2s||zpH0>i 
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Proof.  By  Schmidt's  inequality  we  have 

/+1  4  r+l 

i  f'<x)2d*  <-2fi i-f  f 2(x)dx 


(3.31) 


when  f(x)  is  a  polynomial  of  degree  not  higher  than  N.  (See  [6]), 
(3.31)  then  yields  (3.30)  easily. 

The  next  lemma  follows  easily  from  the  previous  lemma. 

Lemma  3.3.  Let  z^(x)  £  Pp(Q(l)).  Then  for  any  integral  k  >  0 


'k.Q(l)  -  C(k)p 


2k, 


V'O.QU) 


(3.32) 


Proof.  For  every  fixed  x2  we  have  zp (Xp  x2)£  Pp (I)  and  therefore 
using  Lemma  3.2  we  get 


+1 

2  4  2 

(x1,x2)]  dx1  <  Cp J  zp(x1,x2)dx1 


(3.33) 


Integrating  (3.33)  with  respect  to  x2  we  obtain 

0,Q(1)  -  Cp  I  1  Zp I  1 0,Q(1)  (3>34) 

3z 

and  analogously  for  -^E-.  By  Induction  we  get  (3.32).  Obviously  (3.32) 
is  equivalent  to  (3.29)  and  therefore  Theorem  3.4  is  completely  proven. 

3.5  The  convergence  of  the  p-version  of  the  finite  element  method 

Theorems  3.1  and  3.2  lead  immediately  to  an  a  priori  estimate  of  the 


rate  of  convergence  of  the  p-version  of  the  finite  element  method. 


THEOREM  3.5.  Let  £  H  (^g),  k  >  1  be  the  exact  solution  of  the  problem 
(3.1),  (3.2)  and  let  u^  be  the  finite  element  approximation  then 


K"up^i,n0  5  c<k’t)p"<k~1>+£IKHk,n 


(3.35) 


when  £  >  0  is  arbitrary.  For  the  boundary  condition  T  =  ,  £  can  be 

u  dn 

set  equal  to  zero. 

2 

A  polynomial  of  degree  p  has  N  degrees  of  freedom  with  N  Ps  p 

therefore  (and  P^I)  has  dimensi'on  N  with  N  p^  and  (3.35)  can  be 

P  P  >  U 

rewritten  in  the  form 


(k-l) 


+e 


I  VUP*  'l.fjg  -  C(k’  £)N 


o' 'k,n0  • 


(3.36) 


For  the  conventional  finite  element  (h-version)  approach  with 
quasi-uniform  mesh  we  have 


'  VV  li,a0  i  c  h  lluoMk,o0 


(3.37) 


with  p  =  min  (k-1,  q) 


where  q  is  the  degree  of  the  complete  polynomial  used  in  the  elements. 

Realizing  that  in  this  case  the  number  of  degrees  of  freedom  N  satisfies 
-2 

N  »h  we  can  rewrite  (3.3.7)  in  another  form 
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and  this  rate  of  convergence  is  an  optimal  one  (possibly  up  to  e  >  0 
arbitrary)  (see  [3]).  This  shows  that  the  p-version  gives  results  which 
are  (neglecting  £  >  0)  not  worse  than  the  conventional  h-version  with 
quasi-uniform  mesh  if  we  compare  the  number  of  degrees  of  freedom  leading 
to  the  same  accuracy.  In  addition  the  convergence  can  be  much  better 
because  we  do  not  have  the  restriction  on  the  convergence  rate  due  to  the 
degree  of  the  elements  as  we  have  in  the  usual  h-version. 

Further  as  will  be  seen  in  the  next  section  (see  theorem  4.3)  under 
some  conditions  which  are  usually  satdsfied  in  practice  the  factor  1/2  in 
(3.36)  can  be  removed  and  then  the  p-version  will  be  superior  in 
comparison  to  the  usual  (h-version)  finite  element  method  with  quasi¬ 
uniform  mesh. 

Let  us  remark  on  the  other  hand  that  when  the  usual  (h-version) 
with  the  proper  refinement  of  elements  is  used  then  in  general  the 
convergence  rate  can  be  better  than  in  the  case  of  the  p-version  with 
fixed  mesh  -  see  [3].  Although  the  general  theory  for  a  method  combining 
the  h  and  p  version  in  an  obvious  manner  is  not  yet  developed,  we  can 
expect  that  the  theoretical  and  practical  advantages  of  both  approaches 
can  be  combined. 

Let  us  assume  now  that  the  convergence  rate  of  the  p-version  of  the 
finite  element  method  is  r,  i.e.  assume  that 


u  -u 


i,n 


<  K 


-r 

P 


(3.39) 


Then  the  following  theorem  holds. 
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THEOREM  3.6.  Let  Uq £  H^ft)  and  assume  that  (3.39)  holds.  Then 

1+r-e  — 

i)  Ug  G  H  (ft*)  where  0*  is  any  domain  such  that  0*  C  T^, 

i  *  1,...  m  where  T  are  the  triangles  of  the  triangulation  S  and 

I  l"0l  I !«-€,<!*  i  Mfl..r.e)(||»0||1>a()+O  (3.40) 

ii)  uQe  Hi+r^~e(Tt),  1  -  1, ...» 
and 

^u0^1+r/2-e,Ti  -  A(Ti>r>£)(|  lu0l  li,flo+K^  (3.41) 

Proof.  Theorems  3.3  and  3.4  are  obviously  valid  not  only  for  a 
rectangle  Q  but  for  any  parallelogram. 

i)  From  theorem  3.3  we  see  that  (3.42)  holds  for  any  ft*  of  the 
form  of  a  parallelogram.  This  is  obviously  sufficient  for  (3.40)  in 
general. 

ii)  Because  any  T^  can  be  covered  (with  overlapping)  by  a  finite 
set  of  parallelograms  (3.41)  follows  directly  from  Theorem  3.4. 

The  practical  importance  of  Theorem  3.6  lies  in  the  observation 
that  the  triangulation  of  ft  has  to  be  made  so  that  the  possible  singu¬ 
larities  are  located  at  the  boundaries  of  T^ .  Exactly  this  was  done  in 
the  linear  elastic  fracture  mechanics  problems  analyzed  by  Szabo  and 
Mehta  [26]  using  the  p-version  of  the  finite  element  method. 
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4.  THE  SINGULARITY  PROBLEM  AND  THE  P-VERSION  OF  THE  FINITE  ELEMENT  METHOD 

4.1  Preliminaries 

In  this  section  we  will  write  §  instead  of  Q(l).  §(p)  was  defined  in 
section  2  .  Let  T^  be  an  open  triangle  with  the  vertex  in  the  origin  and 
T  C Q  (1/3)U(0,0).  See  figure  4.1. 


4.1  Triangle  with  vertex  at  singularity 

Denote  the  sides  of  T  going  through  the  origin  by  s^  and  s2  and  the  remaining 
one  by  s^- 

A  A  A 

By  $  we  denote  the  mapping  (one  to  one)  of  Q(ir/2)  onto  Q,  defined  so  that 

HO  *  x,  5  i  (5r  52)ak*/2) 

x  =  (x^,  x2)EQ 
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with 


•  2  r 

xi  -  sin  C± 


1,2. 


By  4  we  denote  the  inverse  mapping  to  4.  Further  let 


<i>  ^—1 

T  =  4  (T), 


4  ~  -i 

s .  =  4  1  (s  . ) 
i  i 


i  =  1,2,3. 


.4 


T  now  will  be  a  curvilinear  triangle  with  smooth  sides  and  positive 

2  2 

angles.  In  fact  the  line  x?=  cx.^  (0<c<^)  will  be  mapped  into  sin-£2  =  c  sin  ^ 


and  so 


-  •  1/2  ,  r 

Z, £  =  arcsin  c  sin 


$  A  /  ^ 

Therefore,  T  is  a  curvilinear  triangle  and  T  CQ(pJ  U  (0,0),  o.  »  arcsin  — 

0  0  /T 


We  see  also  readily  that 
sin25. 


sin25- 


is  a  function  bounded  from  above  and  below  on  T 

.1 


4 


Now  let  v  6  H  (T)  be  given  and  define 
V(£)  =  v($(0). 

We  Drove 


Lemma  4.L  Let  v  E  H1(T)  then  V  £  H1(T<t)  and 


Cl  IMI1T  <  II V  ||ltT*  <  c2 


"l,T 


(4.1) 


with  0  <  c^  <  c2  <  00  independent  of  v. 


Proof.  First  let  us  show  that 


* 


-23- 


=1  II  fir"  0.1  i  II  o.T»  id2  '1ST1'  0,1. 


We  have 


Therefore, 


3V  3V  3xi  3V  - 

%  <«>  -  ^rL  a?;  ■  ^  <*«»  si»2V 


/.  <l7>2  «i«2  •  / (I5t>2  s1“2z5i  *k-2  ■ 

T9  ...  T 

sin25 

Because  as  we  mentioned  — : — rr-  is  a  function  bounded  from  above  and  below,  we 


sin2£„ 


get  (4.2)  , 
Further  we  have, 


/■  2  ^1 ^2 _ 

JT  V  sin2CL  sin2C2 

[jfv  p<ix] p  [-4  (sin^i 


-)q  dx 


i  +  -  =  1. 
p  q 


Because  in  the  neighborhood  of  the  origin  we  have  it  follows  that 


- i —  -a*  x-1^2  and  therefore  for  p=3  and  q=3/2,  the  second  term  in  (4.3)  is 

sin2£^  i 

bounded-  On  the  other  hand  by  the  Sobolev  imbedding  theorem  we  have 
f  ,  1/3  2 


and  so  we  get 


dx  <  C  v 


«  vd«  <c||v||ltT< 


Now  (4.4)  together  with  (4.2)  gives 
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(4.3)  yields  also 

INI  0,1  ill'Uo,!4 

and  we  easily  complete  the  proof  of  the  inequality  (4.1) 

A 

Lemma  4.2.  Let  v(x)  be  defined  on  1(1)  (as  a  function  of  one  variable 
0<x<l)  and  let 


/  v2 ,1~1  dx  +  l  ‘S’ 


2  2 
x  dx  <  A  <  00 . 


Let  S  be  the  triangle  with  vertices  (0,0),  (1,0),  (1,1)  (as  in  Figure  3.3) 
let 


Then 


u(xL,x2)  =  v(x1)(l  -  — ) 


u'l  i,s  1CA 


with  C  independent  of  v. 
Proof.  We  have 


fa  “2te  '  /  v2<xi)‘,,‘i  _/*(1  -  ^>2 

<_  2J  v2(x1)  dxl  [x^  +  y  x^]  £ 


CA 


.(4.5) 

and 


(4.6) 


(4.7) 
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Fur  t  her 


f^2  *  <-  2[/^2  -i|4'2-2  +P  7^  ** 


-  2[J  +/h'V 


<  CA  (4.8) 


/(l^)2dx  1  ^  v2dxl  /  - 

S  1  0  X1 

i  /v2  x^1  dx1  <  A2 


(4.9) 


Combining  (4.7)  (4.8)  (4.9)  we  gee  (4.6)  and  the  lemma  is  proven. 

Remark.  If  v  is  a  smooth  function  (4.5)  implies  that  v(0)  =  0.  In  addition 

if  v  is  a  polynomial  then  u  is  a  polynomial  (in  two  variables)  too. 

4.2  Approximation  Properties  of  the  space  P  (T) 

P 

Let  us  introduce  a  one  parameter  family  ^(A)  0<A<A^  (y>0,  fixed)  of 

functions  defined  on  Q.  A  function  uA(x)£  V  (A)  iff 
/\  *• 

i)  u^  G  E(Q)  (not  E(Q) ) 

ii)  Supp  u^R^,  a  >  1 

where  we  define 

/V  x.. 

Ra  i  (xGQ(l/3),  <x2<ctx1} 

iii)  It^uJ  <_  C(|k| )  ]  x  [ 

for  ]  x  (  >  A  and  any  k  =  (k^,,^),  ^  _>  0  integral  with 
]  x  [  »  min  (x^  x2) 


{a} 


a  for  a  ^  0 
0  for  a  <  0 
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iv)  u^|  <_  C(|k|)  A  fi’  an£j  q (k)  j_s  independent  of  A. 

Denote  now  as  before 

UA  (£)  =  uA  ($(£)) 

and 

^  (A)  =  {UA(0  |  uA  E  4-y(A)}. 

Mow  we  prove  the  following  theorem: 

THEOREM  4.1 
Let 

UA(0  E  f*  (A) 

Then 


i 

I 

I 

I 

I 


A 


and 


UA  £  H  (Q(ir/2))  for  any  k>0 

-{l/2{k-2y}-l/2} 


H  V  k.Qdr/2)  1  C(k)A 

with  C  independent  of  A. 

First  we  introduce  some  auxiliary  lemmas. 

Lemma  4.3.  For  0<t<7r/2,  and  n  1,  l<k<n,  k,n  integral  define 


\(t)  =2  (-l)k'j(l>  sin2(k'j)t-^  sin2jt 


then 


3-1 


|\  (t)  1  <  C(n)  t{2k"n} 


dt 


(4.10) 


(4.11) 


(4.12) 


Proof.  Obviously  (t)  is  a  trigonometric  polynomial.  In  the  neighborhood 


of  t=0  we  have 


,  2k  2k  ^  .  ,  2k+2, 
sin  t  =  t  +  0  (t  ) 


’’X 
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and  therefore. 


,n  .  2k  r , 

-  (t)  |  <  C(n)t{2k'n} 

dt 


Hence  for  j  <  k  <  n 


(4.13) 


si„2(k-i>  t  <  c(n)c2(k_:)>  t(2i->> 

dtn 

<  C (n)  tAi(k,:i,n) 


where 


A1(k,j,n)  =  2 (k-j )  +  {2j-n}/ 

It  is  easy  to  check  that 

A^(k,  j  ,n)  _>  {2k-n} 

and  this  yields  (4.12). 

Lemma .  4.4 
Let 


Then 


uA  e  r  (A) 

3  !  k  1  u 

I — ~k—  \  (O I  <  c(|k|)  ]?[  _{|k|~2Y} 
8C1lcl3C2  2 

1/2 


for  1  T  [  >  6,  0=arcsin  A 


(4.14) 


(4.15) 


Proof.  We  have  (see  e.g.  [9],  p.  19) 


3klUA  T^k  ,  9Jua 


(*(5)) 


i 

I 


(A. 20) 


By  simple  computation  we  get 


A2  >  -{ | k| -2y} 

and  now  (4.15)  follows  from  (4.19)  and  (4.21). 

Now  we  will  prove  Theorem  4.1. 

Proof  of  Theorem  4.1.  Define 

$  "-l 

a  1  a 

$ 

3ecause  by  assumption  supp  D^C^,  we  have 

•'  UA'Lq(tt/2)  -  '  UA'k,Q(e)  Or$  +  H  Ua|£,  R4 

a  a 


Now  we  will  separately  estimate  the  terms  in  (4.22). 
On  Q(Q)  using  (4.16),  lemma  4.3, and  the  property  (iv)  of  u^ 

'ml  i 

1  'U  1  2 

2  U.A  I  ‘ 


- — -  <  C(m  ,  m.) 

3e1"1352”2l  "  2-  2 


mi  n>2  2{2j-m, }  2{?X-raJ 

-  2E  h  1  ’2  2 

j=i  £=i 


1/2 

Because  5  <_  CA  we  get 


-  Q(3) 


-2{j+£- 

A 
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3'm|U, 


m2 


0*Q (@  ) 


mi  m2 


<  C(m1,m2)  22  A 

j=l£=l 


-2{j+  l- y}  +  1/2 [ 2{2 j-m  }  +  2{2-£-m2}]  +  1 


m,  m. 


<  C(m1,m2)22 

j-U-l 


9  A_  (j  ,£,m  .m^y)  +  1 
A  1  l 


<_  C  (tn^ ,  m2 )  A 


-{m-2y}+l  ' 


(4.23) 


On  -  Q (9)  using  lemma  4.4,  ve  have 
9  I  ml  U.  1,2 


■o,R2  -  Q(0) 


<  C  (m 


it/2 


rm2)  f  h  h 

Jl/2 

CA 


-2{ I  ml -2y} 


<  C(.lt  .2)  ca1/25r{2tl"l-2v>-2>. 


(4.24) 


(4.23)  and  (4.24)  yield  (4.10)  . 


THEOREM  4.2.  Let  6  f  (A),  be  continuous  on  Q,  u.(0,0)  =  0,  and  u  =  0  on  the  side 


s^  of  T.  Then  there  exists  z  €  P  (T)  such  that  for  any  k  2y+l,  k  integral 

P  ■&- 


i) 


uA~zp"l,T  —  C(k)p 


-2)  .-1/2  (k-2y}-l/2 


(4.25) 


ii)  z  =  0  on  the  side  s.  of  T 
P  J 


n 


»  .  j 
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/  <V  V'  ai_1  *  j  &  (^p))2ll  ^ 

i  s .  1 

L  i  J 


<_  C(k)p-^k-2^  £~l-/2{k-2y}+l/2'  ^,2  <4-26> 

where  we  denote  by  the  length  parameter  of  the  side  s.^  measured  from  the  origin. 

Proof.  By  theorem  4.1  the  function  satisfies  (4.10).  Therefore  there  exists 
a  sequence  of  trigonometrical  polynomials  with  period  it  and  symmetric  with  respect 
to  the  lines  ^  =  ±  tt/2,  ^  =  0,  =  ±  it/2,  ^  =  0  such  that  for  0  <  m  <  k 


V'pH.,  Q(,/2)  i  C(".WP'<k'">  d'l/2{1<'2Y,+  1/2' 


(4.27) 


Recalling  that  U.  =  0  at  the  vertices  of  T  we  can  modify  t  so  that  t  =  0  at 
A  P  P 

these  vertices  and  (4. 27)  holds  for  all  m  >  1  +  e.  In  addition  using  the  trace 
theorem  we  have 


Wll,  s*  £  C(k)p"(k-2)  A'1/2{k"2Y}+1/2  i=l,2,3. 


(4.28) 


Defining  zp  ($(£))  *  tp(£)>  zp  is  an  algebraic  polynomial  of  degree  p.  Using  lemma 
4.1  we  get  (4.25).  By  assumption  tp  =  0  and  =  0  at  the  vertices  of  T$  and  there¬ 
fore  zp  =  0  at  the  vertices  of  T.  Further  we  have 

i  1  ‘  c  (4?>1/2  11  V"»“i,  ,J 

where  we  have  denoted  the  length  parameter  of  s?,  i  =  1,  2  measured  from  the  origin 
$ 

by  -4  .  Therefore  for  z  (4.)  on  s.  we  get  for  i  =  1,2 
i  P  l  i 
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J  (z  -u)2  4  .~1  di  .  <  C  f  (t  -D,)2(4  V2  ^  d 

r  P  1  1  ~  A  P  A  1  1  1 


1°  f 

J  $ 


W2^"1  *! 


<  C  t  -U.  L  <5 
—  "  p  A"l,  s .  ■ 

i 


(A. 29) 


Further, 


/  'a^V^*2  h  “i  i  c  f„  (VV,2iN! 


1 

1 

vvf  s$ 

F  ’  i 


(4.30) 


and  therefore  canbining  (4 . 29),(4 . 30)  and  (4.28)  we  get  (4.26). 

A 

Realizing  that  on  the  mapping  $  is  an  analytic  one  we  can  use  lemma  (3.1)  and 
achieve  z^  =  0  on  s^.  So  theorem  4.2  is  completely  proven. 

Remark  to  theorem  4.2,  We  have  assumed  that  triangle  T  is  situated 
as  in  Fig.  4.1.  It  is  easy  to  see  that  a  linear  transformation  of  the  coordinates 
does  not  make  any  change  in  the  theorem.  Therefore,  theorem  4.2  is  true  for 
any  triangle  T  with  vertex  at  the  origin. 

4.3.  A  concrete  family  ^(A). 

First  denote  by  X (x)  0<x<°°  a  function  with  all  continuous  derivatives 

such  that  X (x)  =  0  for  0  <  x  <  1/2,  and  X (x)  =  1  for  l<x<°°.  Further  let 
XA(x)  -  X(|). 
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Let  now  u  £  Q 

u  =  p (r)  0(<j>)  (A. 31) 

be  given  with  r,  <j)  being  the  polar  coordinates  and  0(4))  be  a  function  with 
all  continuous  derivatives.  Further  assume  that  u  has  support  in  and  u  =  0 
on  Sy  In  addition  let  p(r)  be  continuous  P(0)  =  0  and 


^-£>1  <  rY~n  c(n) 


dr 


n'  — 


with  y  >  0 


(A. 32) 


Now  let  P^=  X^(r)p(r)  and  u^=P^0(<Ji).  Then  obviously  u^=  0  on  s^  and  has  compact 
support  in  for  all  A. 

Let  us  show  now  that  u^  is  a  (A)  family  of  functions.  Obviously,  parts 
i)  and  ii  )  of  definition  of  ^(A)  are  satisfied. 


Further 


a  u, 


3r 


t3J  XA  -^-i 
j — 4i  li — p 


ill  lr. 

j=0  3rJ  3r 


_j!  1  A"j  AY~k+:i<c(k)A 

j-0 


Y-k 


(A. 33) 


because  of  (A. 32)  and  the  fact  that  xA  =  0  for  r  1  A/2, 


and 


3 1  k'u 


3x1kl3x2kz 


3£XaP 


t-Q  3r 


IkR 


3jXA  3*~jP 


J"' 


£=0j=0 

W 

1  C  t?~l  A"lkl+£  <  C(|k|)AY 
1=0 


-  k 


(A . 3A) 


so  property  iv)  of  the  'F^(A)  family  is  satisfied. 
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Because  X^(r)  “  1  ^or  r  _>  A,  (4.31)  yields  the  property  iii)  of  the  f  (A) 
family. 

Let  us  now  show  another  property  of  u^  in  our  concrete  case. 

Lemma  4.5.  Let  u  be  given  by  (4.31)  with  p  and  0  satisfying  the  conditions  spelled 


out  above  and  let  u^  =  xa  (r)p(r)Q(<j>) .  Then 


™A.  5  -  c  a 


(4.35) 


with  C  independent  of  A. 

Proof.  We  have  v  =  u  -  u^  =  (l-x^)p(r)0($)  and  therefore  we  have  using  (4.32) 


ifri  i  'Y_1  « 

i 

and  v  *  0  for  r  >  A.  Therefore 


(4.36) 


f  (— )2  dx  <  C 

f  r2y~2  r  dr  <  C  A2Y 

(4.37) 

'  3x 

i 

Jo 

We  have  also 


/ |v|2dx£C ^  r2Y+1dr  £  C  A2y+2 


(4.38) 


Combining  (4.37)  and  (4.38)  we  get  (4.35). 

From  the  point  of  view  of  applications  the  function  p(r)=r  ^  g(|£gr|)  is  of  impor¬ 
tance  especially  with  g(x)  =  x*5  or  g(x)  =  cos  x  etc.  Then  (4.32)  is  satisfied  with 
Y  *  x  -  e,  e  >  0  arbitrary. 

4.4  The  Convergence  Rate  of  the  p-Version  of  the  Finite  Element  Method 

Returning  to  our  model  problem  (3.1)  (3.2)  we  can  assume  (see  e.g  .  [10],  [15]  that 
its  solution  u^  can  be  written  in  the  form 


u0  =  0)  + 


LX 

r  ^ 


(4.39) 


with  w  £  Hk(ftg)  ^  Hq  (0q)  for  the  boundary  condition  fu  =  u  and  116^  (ftQ)  for 


the  boundary  condition  Tu  *  and 


vi  =  aipi(ri)  9i  (V  e  H0  (V 


(4.40) 


(resp.  H  (^g))  where  are  the  polar  coordinates  with  respect  to  the  vertices  A^. 


of  the  polygon  and  a  are  constants  with 


Pt  (r)  =  rYi  g± ( |  -CgrJ) 


(4.41) 


Ai'" 


1^-!  <  Cl(j)  x  +  V  0  <  x  <  -  PM  >  0,  s  -  1,2..  ,  . 

and  0^  is  a  function  with  all  continuous  and  bounded  derivatives.  The  coefficient 

is  closely  related  to  the  angle  in  the  boundary  of  (2  at  the  vertex  A^.  Without 

any  loss  of  a  generality  we  can  assume  that  0^  are  smooth  periodic  functions  with 

2 

period  2tt so  that  the  function  v^  is  defined  in  the  entire  R  .  This  form  occurs  in 
all  elliptic  problems  e.g.,  elasticity,  see  [15]. 

Let  now  S  be  a  triangularization  of  such  that  all  vertices  of  are 
vertices  of  the  triangulation.  Obviously,  we  can  assume  that  the  support  of 

t 

p^(r)  is  arbitrarily  small  and  v^  has  support  in  an  open  cone  (with  angle  <  tt) 
so  that  the  triangle  T^  lies  inside  such  a  cone.  See  Fig.  4.2. 


Cone  Ki 


<>•  Cone  K2 


VA-LVlXi,A-Z  ”i,A 
i=l  i-1 

where  x*  »  is  Che  function  introduced  in  section  4.3  with  respect  to  the  origin 
(denoted  by  index  i)  of  the  polar  coordinates  of  v^.  Now  we  are  able  to  prove  the 
major  theorem  of  this  section. 

THEOREM  4.3.  Let  u^  be  the  exact  solution  of  the  problem  (3.1)  and  (3.2)  which 

can  be  written  in  the  form  of  (4.39),  (4.40),  (4.41)  and  let  u  be  the  finite 

P 

element  approximation.  Then  for  k  >  1 

II  u0-up  II  1,$}  -  pW  +  £|lfHka,  e>0  arbitrary  (4.42) 

M  *  min  [k-1,  2y^]  (4.43) 

where  y  ”  ~~  and  a  is  the  angle  of  at  the  vertex  A  . 
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Proof.  The  exact  solution  uQ  can  be  written  in  the  form  (4.39)  with 

77 

u  S  H  (Q)  and  y.  *  —  .  See  [15].  It  is  sufficient  to  show  that  functions 
i  a . 

■L  i  i 

cj  and  v.  can  be  approximated  by  z  £  P  _L  J  resp.  z  £  P  1  (for  boundary 
i  rr  P  P,0  P  P 

conditions  Tu  =  u  respectively  Tu  =  -j^)  preserving  the  estimate  (4.42). 

dn 

Using  theorems  3.2  and  3.1  we  see  that  the  function  t)  is  approximable  in  the 


desired  way  and  we  have  to  concentrate  only  on  approximation  of  the  functions 
E  Hq(S)q)  resp.  H^(Sig).  It  is  easy  to  see  that  v.  ^  £  Hg(Qg)  ^ 

v  E  Hg(ftg)  and  using  lemma  4.5  we  see  that 

'K.a  -v11I  1,  n0  i  “•44) 


In  addition  using  theorem  4.2  and  the  remark  to  it,,  there  exist  polynomials 

z  •  E  P  (T  ) 
p,j  7  y 

such  that  for  any  k.  >  2y.  +  1 
i  —  i 

II  zp,j"Vi,All  l,T  1  C(ki)p"(ki'2)  A  "1/2'ki‘2Yi}  +  V2  (4.45) 

and  z  .  satisfies  also  the  condition  (4.26)  on  the  sides  of  T . .  The 
P>  J  J 

polynomials  z  ,  are  not  in  general  continuous  through  the  sides  of  the  triangles 
P  >  J 

T.,  nevertheless  because  of  condition  (4.26)  the  function  z  .  -  z  .  ,  - . 

J  p.jj^  pJ2  deflned  on 

the  common  sides  of  T  and  T.  satisfies  the  condition  (4.26)  too,  it  is 

J1  J2 

a  polynomial  of  degree  p,  and  is  zero  at  the  end  points  of  the  side  s.  Using 

lemma  4.2  we  can  add  a  polynomial  5  .  of  degree  p  on  T.  so  that  the  continuity 

P.j 

through  the  side  s  is  accomplished  and  preserving  the  estimate 


5  -v  II  <  c(k  )p-(ki-2)  -l/2{k  -2y  }  +  1/2 

zPl  ''i.Al'l.Qg  -  C(VP  1  A  11 


(4.46) 


In  addition  if  v,  A  *  0  on  a  side  aC3fL,  then  z  *  0  on  this  side  also  and  so  if 

i ,  A  0  p 

vi e  h5(v  then  vi,AeHJ(v  and  zP.G  PJS1  (v  • 


0 


Combining  now  (4.44)  and  (4.46)  we  see  that 


-V  II  <  C(e>  4V%c(Vp-^-2)a-1/2(V^)  +  1/2 

i  1  l,uo 


Select  now 


ki  ■  V  k0  1 i 


such  that 


and 


-2  +  k„ 


(l/2)k0-(l/2)  - 


>  2-e 


A  =  p 


-X 


Then  we  have 


~2y.  +  e  ^  .  -kn  +  2  -y,A+  A{-?  -1/2} 


?±  -^H^l^Pi  i  +C(VP 

<  C(e)P  ~2yi  +  e>  +  c(kn)P~2yi  +  e" 


<  c(e)P'2Yi  +  £ 


(4.47) 


and  so  (4.47)  yields  immediately  (4 . 42) (4 . 43) . 

So  far  we  have  analyzed  the  rate  of  convergence  for  the  model  problem  (3.1) 
(3.2).  It  is  obvious  that  the  model  problem  was  not  an  essential  one.  It  was 
essential  only  to  analyze  the  approximation  behavior  in  H2(Oq)  resp. 

»5  «v  • 

Combining  the  main  result  of  the  theorem  4.3  and  Theorem  3.6  we  see  that 
(up  to  e)  the  estimate  (4. 42), (4. 43)  is  a  best  possible. 
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5.  NUMERICAL  EXAMPLES 

In  order  to  illustrate  the  results  of  the  theorems  and  in  order  to  show 
the  efficiency  of  p-version  of  the  finite  element,  we  now  present  several 
examples.  The  first  example  is  a  simple  bar  problem  in  one  dimension  and  the 
numerical  results  are  based  on  a  computer  program  written  specifically  for  this 
problem.  The  other  examples  are  two  dimensional  and  the  numerical  results  are 
based  on  COMET-X,  an  experimental  prototype  for  a  general  purpose  finite  ele¬ 
ment  computer  program  developed  at  Washington  University  which  implements  the 
p-version  of  the  finite  element  method  {,  2  ] . 

5.1.  A  One  Dimensional  (Bar)  Problem 

We  consider  the  problem:  ’  *  ^  ,  S3  *  (-1,1) 

u"  -  -q(x)  for  x  e  S3  (5.1) 
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In  this  one  dimensional  case  it  is  possible  to  prove  the  direct  and 
inverse  approximation  theorems  by  using  the  weighted  Sobolev  (respectively 
Besov)  spaces  associated  with  the  Legendre  differential  equation: 
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If  we  let  U  = 
in  strain  energy. 


u| |  denote  the  strain  energy  then  U-U  = 

t  p 


P  E 


is  the  error 


Case  A  du  _  J  2  .  d  2 

to  _Vl"X  *  q(x)  "  "  toVl_X 


1  /  2  -1 

In  this  case  u(x)  =  —  (x\l-x  +  sin 

u(-l)  =  -  f,  u(l)  =  j  . 


x)  and  the  boundary  conditions  are 


Also  the  energy  is 


(l-x~)dx  =  -j  . 


The  coefficients  a^  in  (5.5)  can  be  evaluated  explicitly.  First  (5.5)  becomes 
2i+l  r1  J~  2  „  .  ,  J 

a±  =  — —  J  Vl-x  PL(x)dx  (5.8) 


How  a^  =  0  for  i  odd,  and  using  the  recurrence  relation  for  derivatives  of 
Legendre  polynomials  [1] 

Pi+l(x)  “  PI^l(x)  =  (2i+l)P.(x) 
for  i=2m  m=l,2,...  we  obtain 


f< 


1-x  P2m(x)dx 


1 

_ L_  f 

4nH-l  y. 


11  VL7 


(P2nrfl(x) 


P2m-l(x))dx 


i  r 

JSiJ  co! 


cos  9  (P2aH-l(cos  9)_P'>~.1  (cos  9^)d9  <5*9) 


2m-l 


From  [1],  page  785  formula  (22.13.7)  we  have 

I.  <c°s  9)fwc°»  ^  -  -251  (2:)(2£i )  • 


(5.10) 
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Substituting  (5.9)  and  (5.10)  into  (5.8)  we  obtain  through  straightfor¬ 
ward  caluclation 


a 


2m 


2  (2m) +1  _ t _  /  2m  j 

2  42m2(m+l)(2m-l)  ' m  ' 


Using  Stirling's  formula  it  follows  that 


so  that 


2m 


as  m  ■+•  °°. 


(5.11) 


Therefore,  the  square  of  energy  of  the  error  in  (5.7)  is  given  by 


(5.12) 


where  N  denotes  the  number  of  degrees  freedom  (p=N  is  one  dimension) . 

On  the  other  hand,  in  order  to  study  the  convergence  of  the  (usual) 

2 

h-version  with  N  linear  uniformly  distributed  elements,  let  x  =  -1  +  — 

i  N 

i  =*0,1,2, . . .N  and  let  u^(x)  denote  the  corresponding  finite  element  solution. 
Then, 

u  (x  )  **  u(x  )  i=0,l,2,...N 

hi  i 

and  we  can  compute  the  norm  of  the  error  e^  =>u(x)  -  u^(x) 
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We  get  (for  linear  elements) 

For  quadratic  and  higher  elements  we  get 

■>(£) 

(In  this  respect  see  [3]) 

Figure  5.1  shows  in  the  log  scale  the  behavior  of  the  square  of  the  energy 

error.  We  see  that  in  the  case  of  the  p-version  the  rate  is  practically  4 

as  follows  from  the  asymptotic  analysis.  In  the  case  of  the  h-version  the 

asymptotic  range  is  not  achieved  and  we  see  the  rate  about  1.81  instead  2 

2 

Case  B  u(x)  =  |x|3/2  (1-x2) ,  q(x)  =  -  ~  ( |x| 3/2  (1-x2)) 

dx 

The  boundary  conditions  are  u(-l)  =  u(l)  =  0.  The  only  qualitative  difference 
between  this  case  and  case  A  is  that  the  square  root  singularity  in  u' (x) 
now  occurs  in  the  interior  of  Q  instead  of  at  its  boundary. 

We  again  consider  one  interval  using  the  same  basis  functions  as  before. 
(5.5)  now  becomes  , 

ai  ~T~J_1  ZZ  (!*l  (1-x  ))Pi(x)dx 

2 id-1  fl  ,  ,1/2  .3  7  2 

=  2  -1  !xl  ^2~2X  (sign  x)Pi(x)dx 

|  C  if  i  is  even 

)  21+1  fl  1/2  ,,  ,  2.  .  .  . 

I  — 2  JQ  x  (3-7x  )Pi(x)dx 


(5.13) 


(5.14) 


if  i  is  odd 


(5.15) 


of  the  energy  error  vs.  reciprocal  of  the  number  of  degrees  of  freedom 
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Using  Stirling's  formula  it  follows  that  for  i  odd 
a.  =  0(i)  as  i 

Therefore,  the  square  of  the  energy  of  the  error  is  now  given  by 


E 

i=p+l 


2  2 
ai  2i+l 


which  has  the  same  rate  of  convergence  (up  to  log  term)  as  the  square  of  the 
error  lle^llg  for  the  h-version.  This  illustrates  the  importance  of  the  state¬ 
ment  made  at  the  end  of  section  4  that  in  order  to  get  the  full  power  of  the 
p-version,  singularities  must  be  located  at  vertices  of  the  finite  element  mesh. 

To  illustrate  this  point  further,  we  plot  in  figures  5.2,  |  le^l  |2  and 
!  |e  I  lE  for  case  B,  using  one,  two  and  three  equal  intervals  for  the  p-version 

of  the  finite  element  solution.  The  results  are  summarized  in  Table  5.1.  The 

2  1 

convergence  of  the  h-version  remains  the  same  (  |  |e,  | |„  =  0(  —  lg  N)  for  both 

h  E  N2 

cases  A  and  B.  In  case  A  the  convergence  of  the  p-version  remains  the  same 

regardless  of  the  number  of  intervals  (||e||  =  0(-~ ))  whereas  in  Case  B,  the 

E  n2 


order  is  2  for  two  intervals,  whereas  it  is  only  1  for  both  one  and  three 
intervals.  This  is,  of  course,  because  in  case  B  for  two  intervals  the  singu¬ 
larity  is  at  a  vertex  of  the  mesh  whereas  for  one  and  three  intervals  it  is  in 
the  interior  of  elements  of  the  mesh,  with  the  consequent  degrading  of  rate  of 
convergence.  In  case  A  the  singularity  is  always  at  a  vertex  of  a  mesh.  We 
mentioned  here  only  the  case  of  the  h-version  with  uniform  mesh  spacing.  It  can 
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(a)  Singularity  in  the  interior  of  the  element  (b)  Singularity  at  element  boundary 
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TABLE  5.1 


Rates  of  Convergence  for  the  h-Version  (Linear  Elements)  and  the  p-Version  of 
the  Finite  Element  Method  in  a  Bar  Problem 
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Chat  when  an  optimal  nonuniform  mesh  spacing  for  elements  of  degree  p(fixed) 

is  used  then  Me,!  |2  =0  (—(where  function  0  is  independent  of  N,  but  depends 
h  E  \NP/ 

on  p.  In  this  very  special  case,  it  is  possible,  of  course,  to  analyze  in 
more  detail  the  combined  h-p-version,  but  we  shall  not  go  into  that. 

5.2.  Two  Dimensional  Problems  -  An  Edge  Cracked  Panel  and  a  Parabolically 
Loaded  Panel 

We  now  consider  two  problems  taken  from  two-dimensional  linear  elasticity. 
One  is  an  edge  cracked  rectangular  panel',  shown  in  Figure  5.3,  the  other  is  the 
parabolically  loaded  square  panel,  shown  in  Figure  5.4.  In  both  bases  the 
displacement  field  is  of  the  form  u  =  rai}>(0),  where  r  and  9  are  polar  coordinates 
and  $  is  a  smooth  function.  In  the  case  of  the  edge  cracked  panel  a  =  when 
r  is  measured  from  the  crack  tip;  in  the  case  of  the  parabolically  loaded  panel 
a  =  2.74  when  r  is  measured  from  the  comer  of  the  panel  (See  [27].).  The 
computations  were  performed  with  the  computer  program  COMET-X  which  allows  the 
polynomial  order  p  to  be  varied  between  1  and  8.  We  wish  to  illustrate  the 
following  points: 

(a)  As  claimed  by  the  theoretical  results,  the  rate  of  convergence  is 

0  -  U  -  CN~2a  (5.16) 


(when  neglecting  e  and  the  fact  that  the 
domain).  In  Figure  5.5  we  plotted  U  -  U 

P 

edge  cracked  panel  for  two  x/a  ratios.  U 
energy  value  (of  the  half  panel)  obtained 


edge  cracked  panel  is  not  a  Lipschitzian 
vs  N  ^  on  log  -  log  scale  for  the 
is  an  estimate  of  the  exact  strain 
by  extrapolation  from  the  following 


(U-Un)E 


«  '0.6 


8  7  6 


log(1/N) 


FIGURE  5.5 

Edge  cracked  rectangular  panel.  Estimated  error 
in  strain  energy  vs.  reciprocal  of  the  number  of 
degrees  of  freedom 


expression: 


U  =  max 

„  x  , 
0<— <4 
a 


V8  -  U7N7 
N8“N7 


2  2r 

17.385486 

E 


(5.17) 


in  which  the  subscripts  indicate  polynomial  orders,  a  is  the  applied  stress , 
t  is  the  panel  thickness,  E  is  the  modulus  of  elasticity.  Poisson's  ratio 
was  0.3  in  all  computations.  It  is  seen  that  the  slopes  of  the  log(U  -  U^) 
curves  rapidly  approach  2a  =  1.  Significantly,  the  asymptotic  range  is 
entered  at  low,  computable  p  values.  This  has  been  utilized  in  practical 
computations  [26  ] .  A  similar  behaviour  is  observed  for  the  parabolically 
loaded  square  panel  in  Figure  5.6.  Here  the  slope  of  log(U  -  U^)  approaches 
2a  =  5.48.  For  this  problem  a  series  solution  is  available  and  U  can  be 
computed  to  arbitrary  precision  [3  ]. 

(b)  When  the  singularity  is  not  located  at  a  vertex,  the  rate  of  con¬ 
vergence  decreases.  To  illustrate  this  feature,  we  varied  the  parameter  x  for 
the  edge  cracked  panel  (Figure  5.5)  and  computed  2a  in  (5.16)  from  the  7th  and 
8th  order  approximations: 

u-u„ 


log 


8 


2a 


7-8 


log  1 
6  »8 


(5.18) 


X 

for  various  —  ratios.  The  results  are  plotted  in  Figure  5.7.  It  is  seen  that 
2a^_g  decreases  as  the  interelement  boundary  approaches  the  crack  tip  C. 

It  was  found  that  aspect  ratios  as  high  as  300  could  be  employed  without 
encountering  numerical  instability. 


.0 


-1.6 


.2 


-2 


-0.8  -0 


log(l/N) 

FIGURE  5.6 

Parabolically  loaded  square  panel.  Error  in 
strain  energy  vs.  reciprocal  of  the  number  of 
degrees  of  freedom 


FIGURE  5.7 

Edge  cracked  panel.  Variation  of  the  rate  of 
convergence  with  x/a 
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5.3.  Round-off  error 

When  high  order  polynomials  are  used,  the  choice  of  basis  functions 
becomes  important  from  the  point  of  view  of  round-off  error.  It  is  pos¬ 
sible  to  design  stable  basis  functions  on  the  basis  of  theory  developed  mainly 
by  Mikhlin  -  see  [16,  Chapter  2]  and  [4,  Chapter  4,7].  Of  course,  the  choice  of 
basis  functions  is  also  influenced  by  programming  considerations  and  the  range 
of  p  for  which  the  program  is  written.  In  general,  it  is  desirable  that  the 
basis  functions  be  hierarchic,  as  described  in  Section  6.1.,  and  computation  of 
elemental  stiffness  matrices  and  load  vectors  be  as  simple  as  possible. 

The  basis  functions  currently  in  COMET-X  were  chosen  primarily  on  the 
basis  of  programming  considerations  and  they  are  not  optimal  from  the  point 
of  view  of  round-off  error.  Experience  with  the  code  has  not  indicated  signifi¬ 
cant  accumulation  of  round-off  error,  however,  in  double  precision  computations 
within  the  range  of  p  allowed  by  COMET-X  (1  to  8) . 

To  study  the  characteristics  of  these  basis  functions,  from  the  point  of 
view  of  round-off  error,  the  assembly  and  elimination  procedures  were  executed 
in  both  double  and  single  precision  ( 7  resp.  15  decimals  on  the  DEC  System  20 
computer)  for  the  two  problems  described  in  Section  5.2.  All  other  computations 
were  performed  in  double  precision  only.  (COMET-X  employs  a  modified  version 
of  Irons'  frontal  solver  [ll]  to  carry  out  assembly  and  elimination).  The  results, 
given  in  Table  5.2  indicate  that  for  p  <_  8  the  round-off  error  is  not  critical  but 
if  significantly  higher  p  is  to  be  used  then  it  will  be  necessary  to  exercise 
caution  in  selecting  the  basis  functions. 


-57- 


6.  COMPUTER  IMPLEMENTATION  OF  THE  P-VERSION:  COMET-X 

In  order  Co  implement  the  p-version  efficiently  it  is  necessary  to 
have  available  a  family  of  finite  elements  of  arbitrary  polynomial  degree 
having  certain  properties.  The  family  should  allow,  for  example,  as  much 
information  to  be  carried  over  as  possible  when  increasing  the  degree  from 
p  to  p+1.  The  present  version  of  COMET-X  contains  a  family  of  triangular 
finite  elements  which  enforce  C°  continuity  across  interelement  boundaries 
for  problems  which  require  solutions  in  Hq(Q)  (planar  elasticity).  We  now 
describe  some  of  the  salient  feature,s  of  COMET-X. 


6.1  Hierarchic  Property  of  Basis  Functions 

The  basis  functions  corresponding  to  an  approximation  of  degree  p 
constitute  a  subset  of  those  corresponding  to  an  approximation  of  degree  p+1. 
Therefore,  the  stiffness  matrix  of  the  element  of  degree  p  is  embedded  in  the 
stiffness  matrix  of  the  element  of  degree  p+1.  All  calculations  performed 
in  generating  the  pth  order  elemental  stiffness  matrices  and  load  vectors 
can  be  saved  for  use  in  the  (p+l)st  degree  calculation.  We  call  this  the 
hierarchic  property  of  the  family. 

As  an  illustration  of  the  difference  between  conventional  and  hierarchic 
basis  functions,  consider  linear  and  quadratic  C°  basis  functions  for  a 
triangle  (given  in  natural  coordinates  (L^,  L^,  L^) ;  see [20]  for  a  discussion 
of  natural  coordinates) .  The  linear  function  which  is  one  at  vertex  i  and 
zero  at  the  other  two  vertices  is  L  i=l,2,3  and  it  is  the  basis  function 

i 

for  the  nodal  variable  u(i)  i  =  1,2,3.  In  defining  quadratic  approximations, 


1  ■■  v 


conventional  approaches  use  the  nodal  variables  u(i) ,  u(i')  i,i'  =  1,2,3 
where  i'  is  the  midpoint  of  side  i  (opposite  vertex  i) .  It  is  clear  that 
the  basis  functions  corresponding  to  u(i)  i  =  1,2,3  change  from  the  linear 
to  the  quadratic  approximation.  In  the  hierarchic  approach,  the  nodal 
variables  used  for  the  quadratic  approximation  are  u(i),  ugs(i’)  where  the 
subscript  s  denotes  the  dif f erentation  in  the  direction  of  a  side.  For 
p  >  3,  the  external  nodal  variables  used  to  enforce  C°  continuity  are  jth 
order  derivatives  at  the  midpoint  of  each  side  in  the  direction  of  the 
s:.de  3  <  j  <  p.  Other  nodal  variables  (called  internal  nodal  variables)  are 
used  to  complete  the  polynomial  to  one  of  degree  p.  See  [12  ,13,14,19,20,21,22]. 

6.2  Precomputed  Arrays 

It  is  possible  to  compute  certain  elemental  stiffness  submatrices  (cor¬ 
responding  to  a  standard  triangle  [  14  ] )  once  and  for  all,  and  then  to  use  these 
standard  submatrices  in  order  to  calculate  the  element  stiffness  matrices  in 
a  given  problem.  Precomputed  arrays  based  on  hierarchic  families  permit 
convenient  use  of  elements  of  different  polynomial  degrees  in  the  same  mesh 
because  two  elements  of  different  degree  are  easily  matched  along  an  interelement 
boundary.  The  precomputed  standard  submatrices  are  also  hierarchic  in  character  so 
that  one  version  of  these  arrays,  corresponding  to  the  maximum  polynomial 
degree  that  will  be  used,  can  be  easily  stored  on  a  permanent  file.  Precom¬ 
puted  arrays  are  described  in  [23]  and  have  been  incorporated  into  COMET-X. 

6.3  Computational  Cost 

There  are  three  main  phases  in  the  computational  process  of  the  finite 
element  method: 

a)  Input  phase;  which  includes  the  computations  of  elemental  stiffness 
matrices  and  load  vectors; 
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b)  Solution  Phase;  which  comprises  the  assembly  and  elimination 
processes; 

c)  Output  phase,  which  includes  the  computation  of  displacements, 
Stresses,  etc. 

When  the  number  of  degrees  of  freedom  is  progressively  increased,  the 
major  variable  cost  occurs  in  the  solution  phase.  In  a  number  of  numerical 
experiments  performed  with  COMET-X  it  was  found  that  the  CPU  time  for  the 
solution  phase  can  be  closely  approximated  by  an  expression  of  the  form 

g 

a  +  bN  ;  2  <  0  <  2.4,  a  and  b  constants.  Thus,  although  the  stiffness 
matrix  tends  to  be  more  fully  populated  in  the  p-version  than  in  the  h-version, 
sparse  matrix  solution  techniques  have  provided  substantial  reduction  in  the 
number  of  operations  as  compared  with  solvers  which  do  not  account  for 
sparsity  (0=3).  As  has  been  already  noted,  the  solution  technique  in  COMET-X 
is  similar  to  Irons'  frontal  solver  technique. 

Solution  time  information  is  given  in  ELg.  6.1  for  the  edge  cracked 
rectangular  panel  (x/a  =  3) .  The  computations  were  performed  in  double 
precision  on  a  DEC-20  computer,  (DEC  System  2040,  128K  36  bit  word  memory, 
TOPS-20  operating  system) .  The  time  for  the  frontal  solver  includes  both 
the  assembly  and  elimination  procedures.  The  time  is  given  in  both  CPU 
seconds  and  in  Equivalent  Time  Units.  (ETU).  As  in  [23],  an  ETU  is  the 
time  required  for  squaring  a  full  18  x  18  matrix  by  means  of  the  subroutine 
GMPRD  (double  precision)  of  the  IBM  Subroutine  Package.  On  the  DEC-20 
computer  this  operation  requires  approximately  0.33  seconds. 

The  total  time  accounts  for  all  three  phases  of  the  computation,  including 
computation  of  the  displacement  vector  and  stress  tensor  at  six  points  per 


element. 
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6.4  The  h  and  p  versions  of  the  finite  element  method 

Let  us  compare  the  h  and  p  versions  of  the  finite  element  method  on 

the  basis  of  the  present  state  of  theory  and  experience. 

1)  Asymptotic  rate  of  convergence  (in  energy)  with  respect  to  the  number 
of  degrees  of  freedom: 

a)  For  smooth  solutions  the  rate  of  convergence  of  the  p-version  is  not 
limited  by  fixed  polynomial  degree,  as  in  the  h-version. 

b)  In  the  case  of  nonsmooth  solutions,  the  p-version  has  at  least 
the  same  rate  of  convergence  as  the  h-version  (when  the  h-version 
is  based  on  quasiuniform  mesh  refinement)  but  in  practical  cases, 
for  example  when  the  singularity  is  caused  by  corners,  the  rate  of 
convergence  of  the  p-version  is  twice  that  of  the  h-version. 

c)  The  h-version,  coupled  with  optimal  mesh  design,  results  in 
higher  convergence  rate;  however,  the  p-version  can  also  be  used 

in  conjunction  with  optimally  designed  meshes.  In  this  regard,  the 
mesh  design  seems  to  be  much  less  critical  for  the  p-version  than 
for  the  h-version. 

2)  Input:  Because  relatively  few  elements  are  used  in  the  p-version,  the 
volume  of  input  data  is  smaller  for  the  p-version  than  for  the  h-version. 

3)  Round-off:  In  practical  cases  the  round-off  problem  does  not  appear  to 
be  more  critical  for  the  p-version  than  for  the  h-version. 

4)  Flexibility:  From  the  practical,  rather  than  the  theoretical  point  of 
view,  the  flexibility  of  the  p-version  is  somewhat  restricted  by  the  fact 
that  constant  coefficients  are  assumed  over  large  finite  element  domains. 
At  the  present  there  is  insufficient  experience  with  curvilinear  and  other 
numerically  integrated  elements  in  connection  with  the  p-version. 


-62- 


5)  Solution  time:  The  available  experience  indicates  that  for  a  given 
number  of  degrees  of  freedom  the  solution  time  for  the  p-version  is 
about  the  same  as  for  the  h-version. 

6)  Adaptivity:  Development  of  adaptive  finite  element  procedures  has 
now  been  recognized  as  an  important  area  for  research.  (See,  for 
example,  [18J )  .  From  the  point  of  view  of  implementation,  adaptivity 
based  on  the  p-version  appears  to  be  simpler.  Adaptivity  based  on  the 
h-version  poses  difficult  data  management  problems.  See,  for  example 

[5  ,28].  In  principle,  it  is  possible  to  base  adaptivity  on  a  combination 
of  the  h-  and  p-versions  but  such  an  approach  would  again  pose  difficult 
data  management  problems.  A  more  promising  approach  is  to  employ  mesh 
grading  on  a  prior  basis,  either  manually  or  with  standard  mesh  generators, 
and  then  to  make  adaptive  changes  by  means  of  adjusting  p. 
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The  p-version  of  the  finite  element  method  is  a 
new  approach  to  finite  element  analysis  in  which  the 
partition  of  the  domain  is  held  fixed  while  Che  degree 
p  of  approximating  piecewise  polynomials  is  Increased. 
In  this  paper,  two  theorems  are  presented  which  des¬ 
cribe  Che  approximation  properties  of  the  p-version. 

In  particular,  for  the  singularity  problem,  the  p-ver¬ 
sion  has  asymptotically  as  p  -  »  twice  the  order  of 
convergence  of  the  standard  version  of  the  finite  ele¬ 
ment  method,  if  the  number  of  degrees  of  freedom  is 
used  a3  a  measure  of  convergence.  Various  hierarhic 
families  of  finite  elemencs,  designed  for  computation¬ 
ally  efficient  compucer  implementation  of  the  p-ver- 
sian,  are  described.  These  families  include  conform¬ 
ing  C°  and  C1  triangular  families,  and  conforming  C° 
families  or  rectangles  and  tecrahedra. 


Introduction 


In  Che  finite  element  method  the  soluclon  to  a 
certain  type  of  partial  differential  equacion  is 
formulated  as  a  variational  problem.  In  the  conven¬ 
tional  approach. the  solution  is  then  approximated  over 
the  given  domain  by  functions  which  are  piecewise 
polynomials  on  convex  subdomains  (such  as  triangles  or 
rectangles)  and  which  are  globally  in  Ca,  n  ^  0,  where 
n  depends  upon  che  order  of  the  partial  differencial 
equation.  The  degrees  of  Che  approximating  piecewise 
polynomials  are  fixed  (usually  ac  some  low  number  such 
as  2  or  3)  and  the  accuracy  of  che  approximation  is 
increased  by  allowing  h,  che  maximum  diamecer  of  Che 
finice  elemencs,  to  go  to  zero.  We  refer  to  chis 
approach  as  the  h-verslon  of  the  finite  element  method. 
The  h-version  has  been  studied  extensively  and  asymp- 
tocic  error  bounds  js.  h  —  0  are  well  known  for  its 
rate  of  convergence^  - 

In  a  new  appsoach  developed  at  the  Center  of  Com¬ 
putational  Mechanics  at  Washington  University  a  dif¬ 
ferent  point  of  view  is  adopted.  The  given  domain  is 
partitioned  into  convex  subdomains  which  are  kept 
fixed ,  and  che  solution  is  again  approximated  by  func¬ 
tions  -which  are  globally  in  Cn,  n  _>  0,  and  which  are 
polynomials  over  each  convex  subdomain.  Now,  however, 
accuracy  is  Increased  by  allowing  che  degree  p  of  the 
piecewise  polynomials  to  go  to  infinity.  We  call 
this  approach  che  p-version  of  the  finice  element 
method.  The  p-version  is  reminiscent  of  che  classical 
Ricz  method  but  wlch  one  important  difference.  In  che 
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Ritz  method  the  solution  is  approximated  on  the  entire 
domain  by  polynomials  (or  ocher  smooch  functions) 
whereas  in  the  p-version  of  the  finite  element  method 
it  is  approximated  only  over  each  convex  subdomain  by 
a  polynomial  and  globally  che  approximate  is  required 
to  be  in  C“.  This  difference  leads  to  a  rate  of  con¬ 
vergence  for  the  p-version  which  is  higher  than  that 
of  both  the  Ritz  method  and  the  h-version,  and  also  to 
other  computational  advantages. 

Ill,  A  Sample  Problem 

In  order  Co  illustrate  the  application  of  the  p- 
version  of  the  finite  element  method  to  a  practical 
situation,  we  consider  a  sample  problem,  called  Lock¬ 
heed  Test  Problem  No.  2,  which  has  been  used  as  a  test 
case  for  various  finite  element  programs.1-  It  con¬ 
sists  of  a  circular  cylindrical  shell  with  symmetrical¬ 
ly  loaded  cutout3  and  subjected  to  a  uniform  axial  end 
shortening  of  known  amount  (Figure  1) . 


The  3hell  is  made  of  a  homogeneous  isotropic  linearly 
elasclc  material  and  has  constant  thickness,  u,  v,  and 
w  represent  the  longitudinal,  circumferential,  and 
normal  displacements,  respectively.  The  boundary  con¬ 
ditions  at  cne  ends  of  che  shell  are 

w  ■  v  "  ■  0  u  ■  constant  •  0.2  *  10~3  inches 
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finite  element  method.  Then,  ve  describe  several  new 
families  of  finite  elemencs  designed  to  implement  the 
p-version  on  computers.  These  families  have  a  hierar¬ 
chic  property  which  leads  to  compucaclonal  savings 
when  using  che  p-version. 

2.  Theoretical  Background  and  Illustrations 
2.1.  Theoretical  Background 

Lac  ft  be  a  bounded  polygonal  domain  in  two  dimen¬ 
sional  Euclidean  space  R2,  lec  Eift)  be  che  space  of 
all  real  C“  functions  on  ft  with_contlnuous  derivatives 
of  all  orders  on  ft,  let  0(ft)CE(ft)  be  che  subspace  of 
functions  with  compact  support  in  ft,  lac  H°(ft)*L2(ft) 
with  the  inner  product 

(u,v)0,ft  “X  uvdx’  ‘  dxld*2- 

If  ^ 

For  t  >  1  integral  lec  a  (ft)  resp.  Hgift)  be  che  comple¬ 
tions  of  E(ft)  resp.  5(ft)  under  che  norm 


The  10-elamenc  criangulacion  used  co  solve  chi3  problem 
is  shown  in  Figure  2,  together  wlch  che  normal  dis¬ 
placement  along  che  center  line  of  che  shell.  The  num¬ 
bers  4-4-4,  5-5-6,  6-6-7  refer  to  che  degrees  of  che 
polynomials  used  co  approximate  u,v,  and  w.  In  this, 
problem  che  approximations  Co  u  and  v  are  globally  C  , 
and  che  approximation  to  w  is  globally  o'- .  In  Figure  3 
a  comparison  is  given  of  che  numbers  of  elemencs  and 
degrees  of  freedom  used  in  the  application  of  several 
computer  programs  to  chis  problem.  The  name  COMET-CT 
refers  to  an  early  experimental  implementation  of  the 
p-version  of  che  finite  element  mechod.2 
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integral,  1  «  1,2, |a|  ■  a,  +  02-  The  inner  product  in 
H^(ft)  will  be  denoted  by  (*,*)k  g.  For  k  >  0  nonin— 
tegral  Che  space  H^(ft)  and  H§(ftJ*  are  defined  by 
interpolation. 

Consider  the  following  00 cel  problem 
-  iu  +  u  -  f  on  ft,  ft  H°(ft)  (2.1) 


Tu  ■  0  on  3ft 


(2.2) 


where  ft  is  a  bounded  polygonal  domain  and  Tu  *  u  or 
Tu  "■  3u/3n.  We  seek  a  solution  in  che  weak  sense  l.e. 
UqcHq (ft)  such  chac 


B(uq,v)  »  (f,v)  for  all  vcHd(ft) 

(resp.  ueH1(ft))  (2.3) 


where  we  define 


We  remark  that  che  accuracies  of  che  tabulated 
solucions  may  vary  in  any  given  norm.  However,  ic  is 
che  considered  opinion  of  independent  engineering 
groups  chat  the  accuracy  of  each  solution  is  adequace. 
There  is  a  substantial  reduction  in  the  number  of 
finite  elements  which  were  required  when  using  che  p- 
version  of  che  finite  element  method,  to  solve  chis 
problem.  Additional  decails  on  Che  solution  of  Chis 
problem  are  given  by  Rossow  et  al.2 

In  chis  paper  we  consider  che  uniform  (or  quasi- 
uniform)  h-version  and  che  uniform  p-verslon  of  che 
finite  element  mechod  l.e.  che  partition  of  che  domain 
is  refined  uniformly  (or  quasl-unlformly)  while  che 
degree  p  of  approximating  piecewise  polynomials  is 
held  fixed  (in  the  h-version),  whereas  che  partition  is 
held  fixed  while  che  degree  p  is  increased  (in  che  p- 
verslon) .  A  general  theory  for  combining  both  the  h- 
and  p-verslons  in  a  simple  manner  is  noc  yet  developed. 

Early  results  on  che  convergence  of  the  p-version 
of  che  flnice  element  were  empirical  based  largely  on 
numerical  experiments  2. 3. 4, 5, 6, 7.  Recently  8,  a  firm 
mathematical  foundation  was  provided  for  the  p-version, 
in  which  basic  approximation  properties  were  derived. 
We  first  state  two  cheorems  proved  in  8  which  esta¬ 
blish  che  race  of  convergence  of  che  p-verslon  of  che 


a<u0,v)  -  (u0.v)1>fl  (2.4) 

The  concepc  of  convergence  in  che  p-version  of  che 
finite  element  mechod  is  now  formulated  as  follows: 

Lec  S  be  a  (fixed)  triangulation  of  ft,  s  •  (T. } 
id  1,  ._j_.  ,m  where  T^  are  open  criangles  such  chac 
U^i  *  ft  and  T^ ,  Tj  iFJ  have  either  a  coranon  entire 
side  or  a  vertex  or  T^nT.  -  p.  Denote  by 
T'laJCsMa)  the  subset  of  all  functions  ucB^fft)  such 
cnat  if  u/T  ,  is  the  restriction  of  u  to  chen 
u^Xi^ePoCT^T  i.e.  P^Hft)  consists  of  all  functions 
which  are  piecewise  polynomials  of  daggee  ac  most  o  and 
which  belong  to  H1 (ft).  Further,  let 

"  Pj^l  (ft)nHd(ft>.  The  p-version  of  the  finite 
element  method  consists  of  finding  u_  P  ■  1,2,...  where 
UpCpp^g  (ft)  (resp.  P^s‘  (ft))  (for  the  boundary  condi- 
clons'u.-  0  resp.  Tu  ■  (3u/3n)  30  that  (2.3)  holds  for 
all  v*P^  (ft)  (resp.  PlSl  (ft)). 

A  bound  for  che  rate  of  convergence  in  che  p-ver¬ 
sion  of  che  finite  element  is  given  by  the  following, 
theorem: 

Theorems  Lec  ueH*(ft),  k  >  1,  be  the  exact  solu¬ 
tion  of  che  problem  (2.1J,  (2.2)  and  let  u-  be  che 
finice  element  approximation,  chen 


:  “o  -  upii  i.  a  -c(k’e)P 


-Ot-l)  +c  I 


^  k,2  (2'5) 


where  s  >0  is  arbitrary.  For  the  boundary  condition 
7u  *  (3u/3n)»i  can  b«  sec  to  zero. 


freedom  w 
dimension 
the  form 


How,  a  polynomial  of  degree  p  has  H  degrees  of 
;om  with  H  s  p2.  Therefore,  J  and  (?£S(J  )  hi 
ision  .  N  with  N  3  p2  and  (2.5)  can  be  rewritten 


H  u0  "  "pill, a  -COt’£) 


.  (fc-D  +> 


''V'k.ilg  (2.6) 


and  bounded.  Let  $  be  a  criangulacion  of  2  such  that 
all  vertices  of  2  are  vertices  of  the  triangular ion. 

We  now  give  the  rate  of  convergence  of  the  p-version 
of  che  finite  element  method  for  this  situation. 

Theorem:  Suppose  uq,  the  exact  solution  to  the 
problem  (2.1)  (2.2)  can  be  written  in  the  form  (2.9), 
(2.10)  (2.11)  and  let  up  be  the  finite  element  ap¬ 
proximation.  Then,  for  k  >  1 

I  lu0  *  upH  1,0  -  C(e)p  U+  £>0,  arbitrary  (2.12) 


min(k-l,2Yl) 


On  che  ocher  hand,  for  the  h-version  of  the  fin- 
Ice  element  nxechod  vich  quasi-uniform  mesh  we  have 

II  “o  *  “hi!  1,0  -  ch'1l|uollIC>a.  ,J  ■  ■la<k-l*'«> 


where  q  is  the  degree  of  the  approximating  polynomial 
used  in  che  elements.  Recalling  that  in  this  case, 

N  a  h~2  we  can  rewrite  (2.7)  in  the  form 


uo  "  “hH  i,a  -  c  s  2  il  U0H 


This  rate  of  convergence  is  optimal  (possibly  up  to 
t  >0).  Comparing  (2.7)  and  (2.3)  we  see  that  the 
p-version  gives  results  which  are  (neglecting  £  >0) 
aoc  worse  than  che  h-version  with  quasi-uniform  mesh 
if  we  compare  the  number  of  degrees  of  freedom  lead¬ 
ing  to  che  same  race  of  convergence,  Also  the  conver¬ 
gence  rate  can  be  much  becter  because  in  the  p-verslon 
we  do  not  have  che  restriction  due  to  the  degree  of 
the  elements  as  we  have  in  che  h-verslon.  Further, 
in  many  practical  situations,  che  factor  2  in  (2.6) 
can  be  removed,  and  then  the  p-version  will  be  super¬ 
ior  to  che  h-version  with  quasi-uniform  mesh. 

More  specifically  we  consider  the  p-version  when 
used  to  solve  a  singularity  problem.  Assume  that  che 
solution  ug  to  our  model  problem  (2.1),  (2.2)  can  be 
written  in  che  form 

un  *  ,J*  +  Z  vt  (2-9) 

1*1 


with  uiefl  (3)  OHq(3)  for  ch«  boundary  condicion  Tu*u 
and  for  che  boundary  condition  I*u  ■  Ou/3n) , 

’i  -  Vl(V  WC  a>  (reap.  H2(S))  (2.10) 


where  r^,*^  are  polar  coordinates  with  respect  co  che 
vertices  of  the  polygon  2  and  a^  are  constants,  and 


3i(r)  *  ri  «i<U»«f1l> 


|— *  ,J  +D1-  0 


and  3  £  is  a  function  with  all  derivatives  continuous 


It  can  also  be  shown  that  the  estimate  (2.12)  (2.13) 
is  the  best  possible..  Although  the  considerations  which 
lead  co  this  estimate  are  only  for  che  model  problem, 
they  apply  more  generally  co  the  analysis  of  the  be¬ 
havior  of  approximations  co  any  function  in  ai(2) 
resp.  H2(2). 

Thus,  in  terms  of  the  number  of  degrees  of  free¬ 
dom  N,  the  p-version  3olucion  of  the  singularity  pro¬ 
blem  leads  . to  the  estimate 


I  uo  '  V 


s  >  0,  arbitrary 


u  *  min(k-i,2y^) 

whereas  in  the  h-verslon  we  have  the  estimace 


H“0  -  “hll  1.2  i  c  * 


that  is,  the  uniform  p-version  has  twice  che  order  of 
convergence  of  the  uniform  h-version.  Let  us  remark, 
however,  that  when  a  suitable  refinement  of  the 
elements  is  used  in  the  h-version  then  its  convergence 
race  is,  in  general,  becter  than  in  che  case  of  che 
p-version  with  fixed  mesh.  A  general  theory  for  com¬ 
bining  the  theoretical  and  practical  order  advantages 
of  both  che  h-  and  p-versions  is  not  yet  fully  devel¬ 
oped. 

2.2.  Illustrations 

In  order  co  lllustrace  che  results  of  the  theorems 
we  present  cwo  examples.  The  first  is  a  simple  bar 
problem  In  one  dimension  and  che  numerical  resulcs  are 
based  on  a  computer  program  written  specifically  for 
this  problem.  The  second  is  a  problem  In  two  dimension¬ 
al  linear  elasticity  involving  the  analysts  of  a  cen¬ 
trally  cracked  panel.  This  problem  was  solved  using 
COMET- i  (COtiscralnc  METhod-eXperimental)  a  general  pur¬ 
pose  finite  element  computer  program  developed  at  Wash¬ 
ington  University  which  implements  the  p-version  of  che 
finite  elemenc  method. 

2.2.1.  A  One  Dimensional  (bar)  Problem.  Consider 
the  problem:  '  *  ^ ,  2  *  (-1,1) 

u"  *  -q(x)  for  xe2 

where  the  (loading)  function  q(x)  and  the  (Dirichlet) 
boundary  conuudons  will  be  specified  later.  The 
energy  inner  product  is 


B(u,v)  *  (u,v). 


f  u' (x)v' 
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The  weak  solution  ueHq(3)  satisfies 

(u,v)_  "  /  u'  (x)v '  (x)dx  *  /  q(*)v( 

“  <1  •'-l 


in  Case  3 


ii-pii; 


for  all  vtHgCsl) 


First.  we  consider  convergence  when  i 1  is  not  parti¬ 
tioned,  i.e.  ve  use  only  one  interval.  We  choose  as 
basis  functions  1 ,x  and 


V*)  *  pi(c)dt  1  -  1 


where  ?  (t)  is  the  Legendre  polynomial  of  degree  1. 
If  we  write 

u  (x)  *  ■— y—  u(-l)  +■ u(l)  +  V  a.ip.Cx) 

?  2  2  i-1  1  1 


it  follows  from  the  orthogonality  of  Legendre  poly¬ 
nomials  that 

a^  ■  — ^  /  qCxJd^CxJdx  i  -  1,2,...,? 


Also,  denoting  the  error  by  e?(x)  •  u(x)  -  Up(x), 
it  follows  that 

I!  «p!t  E  -  II  »  -  “pll  |  “  II  “II  E  -  ll“p!lg 
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If  we  let  U  •  |  ;u|  |£  denote  the  strain  energy  then 
U  —  tip  —  1 1  1 1  £  is  cne  error  in  strain  energy.  '.Ie 
consider  two  cases: 

Case  A  ^  ,  q(x)  -  -  £  (V77 ), 


u(-l)  -  -u(l)  -  f 
4 


Case  3  u(x)  -  |x|  (1-x2) , 


d2  3/2  2 

q(x)  -  -  — ( |  x |  (l-0),a(-l|  -  u(l)  -  0 


The  qualitative  difference  between  the  cwo  cases  is 
that  in  case  A  the  square  root  singularity  in  u’  is 
at  end  points  of  3,  whereas  in  case  3  ie  is  in  its 
interior.  It  has  been  shownd  that  if  M  denotes  the 
number  of  degrees  of  freedom  (N*p-*-l)  then  as  X  -  » 


This  illustrates  the  importance  of  locating  the  sin¬ 
gularities  at  vertices  of  the  finite  element  mesh  in 
order  to  obtain  the  maximal  rate  of  convergence  in  the 
p-verslon,  as  staced  in  che  second  theorem.  In  order 
to  illustrate  this  point  further  and  to  compare  che 
races  of  convergence  of  che  (uniform)  h-  and  p-versions 
of  the  finite  element  method,  we  plot  in  Figure  4 


ft 


«er 

OwVx*) 


i  ui  w 


j]  ep)|  _  using  one,  two  and,chree  equal  intervals  for 
the  p-version,  and  where  we  have  used  che 

nocaclon  e^(x)  •  u(x)  -  u^(x).  Linear  approximations 
were  used  on  equal,  intervals  in  che  h-version.  It  is 
known  chat  liewllE  J  0((l/}|2)|ln  X|  for  linear  ele¬ 
ments.  and  l!ehiiE  "  0(1/11  for  quadratic  and 
higher  elemencs.  Figure  4  snows  in  the  log  scale 
chat  the  square  of  che  energy  error  in  che  case  of 
che  p-verslon  has  an  exponent  chat  is  practically  4. 
In  che  h-verslon  the  asympcoclc  range  has  not  yec 
been  reached  and  che  race  is  about  1.31  instead  of  2, 
where  3,  che  number  of  degrees  of  freedom,  .is  denoted 
by  NBF,  The  results  are  summarized  in  Figure  3. 
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la  case  a  Che  race  of  convergence  in  che  p-version 
remains  che  3ame  regardless  of  che  number  of  inter¬ 
vals  because  che  singularicy  is  always  at  an  end¬ 
point  of  an  interval,  la  case  3  vnen  the  singularity 
is  at  an  end  point  of  an  interval  (i. a.  when  chere  are 
cvo  intervals)  the  naxiroal  rate  of  convergence  in  che 
p-version  is  achieved,  whereas  when  che  singularity 
is  in  che  interior  of  an  interval  (i.a.  when  these  are 
either  one  or  three  intervals)  the  race  of  convergence 
deteriorates . 

2.1.1.  A  Centrally  Cracked  Panel.  Consider  the 
centrally  cracked  panel  3hovn  in  Figure  6. 


»  >  \  )  »  »  | 
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The  displacements  have  singular  behavior  in  che 
neighborhood  of  the  crack  dp (by  symmetry  only  a 
quarter  of  che  panel  needs  to  analyzed)  in  the  form 
u  *  rfT(.9)  where  r,3  ^are  polar  coordinates  with  re- 
pec  c  to  che  tip,  y  9  7,  and  T  is  a  smooch  function. 
Two  finite  element  triangulations  were  used,  the 
eight  element  mesh  and  che  three  element  mesh  shown 
ia  Figure  6.  In  the  eight  element  mesh  che  polynom¬ 
ial  degrees  were  distributed  In  two  ways:  uniformly 
and  non- uniformly.  In  che  non-uniform  or  graded 
distribution  che  polynomial  degrees  were  greater  chan 
?*3  only  in  crack  dp  elements  (numbered  1,3,4)  and 
che  polynomial  degree  was  held  constant  at  p*3  in  che 
remote  elements  (numbered  5,6,3)  and  che  transition, 
elements  (number  2,7).  For  this  problem  since  y  9  i 


the  theory  predicts  chat  che  error  in  strain  energy 
is 


U  - 


|U“  V'1,1 


"2  -  Cp‘4Y-  C  p'2  -  c  s'1 


if  ve  neglect  the  e  in  che  second  theorem  and  the  fact 
that  che  cracked  panel  is  nod  a  Lipschitzian  domain. 
The  comparable  result  for  the  h- version  is 


»  •  II,,  •  C  S  * 


In  Figure  7  che  computed  strain  energy  0.  (normalized) 
is  plotted  again  X”1  and  che  convergence  pachs  are 
seen  to  be  nearly  linear  for  all  o  >_  i.  Additional 
details  are  given  for  a  cencrally  cracked  panel  by 
Szabo3  and  by  Szabo  and  Mehta12.  In  10  comparative 
plots  of  convergence  in  che  ?-  and  h-versions  are 
presented  for  the  case  of  an  edge— cracked  pant.. 

3.  Hierarchic  Families 

In  order  to  implement  che  p-verslon  of  che  finite 
element  mechod,  it  is  necessary  to  have  available’  a 
family  of  finlce  elements  of  arbitral"/  polynomial 
degree  p.  Although  such  families  of  finice  elements 
have  been  constructed  (by  Kratochvil  ec  ai,11  for 
antample) ,  ve  wish  to  present  a  new  family  vhicn 
the  property  that  when  increasing  che  degree  of  cne 
approximating  polynomial  from  p  co  p-tl  as  much  of  the 
computation  as  possible  is  saved  from  che  pch  degree 
approximation.  This  is  clearly  a  desirable  property 
.or  efficient  computation  when  using  che  p— version  of 
the  finite  element  method.  More  specifically,  in  this 
family,  basis  functions  corresponding  to  an  approxima¬ 
tion  of  degree  p  are  a  subset  of  chose  corresponding 
to  an  approximation  of  degree  jm-1.  Therefore,  the 
stiffness  matrix  of  the  element  of  degree  p  is  a  sub- 
macrix  of  che  sciffness  macrix  of  che  element  of  degree 
p-rl,  and  when  increasing  che  degree  of  approximation 
from  p  co  p+1  only  che  added  rows  and  columns  of  the 
new  sciffness  macrix  have  co  be  computed.  We  call  a 
family  possesing  this  propertv  a  hierarchic  family. 
COMET-X,  che  current  experimental  Implementation  of 
the  p-version  of  the  finite  element  method,  develooed 
at  Washington  University,  is  based  on  cvo  hierarchic 
families  of  conforming  triangular  finice  eleaencs 
for  che  analysis  of  two  dimensional  problems  in  linear 
elasticity.  One  family  enforces  C3  continuity  across 
Interelemenc  boundaries  for  problems  vhich  require 
solutions  In  H r|(.T)  (planar  alascicity) ,  and  the  other 

tan  n«r  is  b*st  iti  nt&M***' 
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enforces  continuity  across  interelement  boundaries 
for  problems  which  require  solutions  in  Kg(^)  (place 
bending) . 

3.1.  Hierarchic  C°  Family  of  Triangular  Elements 

This  family  consists  of  complete  polynomials 
of  degree  p  >  2  defined  either  in  terms  of  nodal 
variables  or” (nodeless)  basis  functions.  The  nodal 
variables  are  divided  into  cwo  classes:  external 
nodal  variables  used  to  enforce  global  C®  continuity 
and  internal  nodal  variables  which  are  added  to  com¬ 
plete  the  polynomial  of  degree  p. 

Let  u  denote  the  approximation  to  the  displace¬ 
ment  field,  and  let  denote  the  deviacive  of  u  in 
direction  of  side  i.  The  external  nodal  variables 
for  each  p  >.  2  are:  u(i),  us  j 'i')  i,  i'  *  1,2,3, 
j  *  2,3,...p  where  i  denotes  vertex  1,  1*  denotes 
Che  midpoint  of  side  i  (opposite  vertex  i) ,  and  us  . 
denotes  the  jth  derivative  of  u  in  direction  of 
side  i. (see  Figure  8). 


"tune  3:  HIBW90CC  C5  TUnratun  E'.adNT,  J  >  2 


For  each  ?  >.  2,  y  (p— 1) (p— 2)  internal  nodal  variables 
defined  as  derivatives  of  order  j,  j  *  3,...,P  evalu¬ 
ated  at  one  vertex,  are  added  to  give  a  complete  poly¬ 
normal  of  degree  p^.  The  basis  functions  for  these 
external  nodal  variables  expressed  in  natural  (tri¬ 
angular)  coordinates  L^,  Lj,  Lj  (see  Feano13,  for 
example,  for  che  definition)  are,  for  vertex  1  and 
side  1 


ca2  -  L3  + 

if  J.  is  even 

L3)Nlj~1)  J  13  °dd 

■mere  (  •  is  the  length  of  side  1.  The  expressions  for 
' '•  pche'r  vertices  and  sides  are  obcaiaed  by  cyclic 

•  ndex  -emulation. 

in  order  to  llluscrace  the  difference  between 
■erarcolc  and  hierarchic  finite  elements,  let  us 
m-ift  --.a  mao  when  ?*2.  For  non-hlerarchic  finite 

•  .eew-i  e  nodal  variables  for  che  linear  approxima- 

"•  .ti.  t  •  1,2,3,  and  for  the  quadratic 
•  .  u'.jn  -he  nodal  variables  u(i')  1  »  1,2,3 
.  -  .  -  «  mu  functions  for  these  nodal  vari- 

.  .«  •  -  e«  t.n  the  linear  approximacion 

•••ponds  to  u(l);  in  the  quadratic 


approximation  L^(2L.  -  l)is  the  basis  function  for 
u(i)  i  *  1,2,3  and  i3  the  basis  function  for 

u(i’)  i  »  1,2,3.  Thus,  the  basis  function  for 
changes  when  increasing  che  degree  of  the  approxima¬ 
ting  polynomial  from  linear  to  quadratic.  For 
hierarchic  elements,  however,  in  the  linear  approxima¬ 
tion  is  the  basis  function  for  u(l)  i  ■  1,2,3  in 
both  che  linear  and  quadratic  approximations  and 
(l/2)LiLi+i  is  the  ba3is  function  for  us2(i’)  in  the 
quadratic  approximacion.  Therefore,  Ehe  basic  func¬ 
tions  for  the  linear  approximation  are  unchanged  wnen 
increasing  che  degree  of  approximation  to  quadratic. 

It  follows  then  that  the  stiffness  matrix  correspond¬ 
ing  to  the  linear  approximation  13  a  submatrix  of  the 
sciffness  matrix  corresponding  to  che  quadratic  approx¬ 
imacion  for  the  hierarchic  family. 

It  Is  also  possible  to  construct  hierarchic  basis 
functions  which  do  not  correspond  to  nodal  variables. 

In  this  case  the  external  nodes  used  to  enforce  global 
C®  continuity  are: 

L^,  Lj,  corresponding  to  vertices  1,2,3 
respectively 

V 

l{L2  -  L1(-L2)j  ,  I^L3  -  4(-L3)J  ,  LJ3L1  -  T.jC-T.j)-!  , 

corresponding  co  sides  1,2,3,  respect¬ 
ively 


1  -  1,2,...,P 


The  internal  basis  functions  are  (j  -2)  independent 
polynomials  each  of  which  contains  factor  Lj^LjCso 
chat  they  vanish  on  che  boundary  of  che  triangle) , 
Cot  J  »  3,...,p  (see  Peano  ^«^4). 

3.2.  Hierarchic  Family  of  Triangular  Elements 


In  this  case  the  situation  is  more  complicated. 

The  hierarchic  triangular  family  described  above 
enforces  global  continuity  (and  no  more  chan  C® 
.continuity,  even  at  vertices).  The  cerm  Constraint 
Method,  in  COMET-X,  is  in  fact  derived  from  the  pro¬ 
perty  that  it  constrains  che  approximating  solution 
co  satisfy  che  degree  of  smoothness  and  no  more  than 
the  degree  of  smoothness  required  by  Che  formulation 
of  the  problem.  In  the  Lockheed  Test  Problem,  (Fig¬ 
ure  1)  an  early  version  of  COMET-X  called  the  Con¬ 
straint  Method  was  used  with  cne  excellent  results 
shown  in  Figures  2  and  3.  It  is  the  property  that  che 
constraint  method  does  not  enforce  more  chan  C*  con¬ 
tinuity  even  at  che  re-entranc  corner  that  contributed 
co  che  good  results.  Although  It  Is  Dosslble  co  con¬ 
struct  a  hierarchic  triangular  family,  it  was  proved 
by  Peano  13,  14 1  chac  in  order  to  enforce  global  C* 
continuity  and  no  more  chan  global  C1  continuity  even 
at  vertices,  certain  additional  constraint  equations 
must  be  satisfied.  Consiu.  for  example,  a  vertex  of 
a  triangular  elemenc  e  and  xec  and  be  coordinates 
along  the  cwo  sides  which  meec  at  che  vertex  (see 
Figure  9) , 
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(1)  A  specially  devised  global  assembly  process 
which  reduces  Che  assembly  of  elements  to  a  stan¬ 
dard  finite  element  assembly  procedure 

(2)  Creation  of  super  elements  (or  macro  ele¬ 
ments)  of  arbitrary  degree  p  >  5.  In  these  macro 
elements  constraints  are  satisfied  within  the 
macro  element  leaving  external  nodal  variables  on 
the  boundary  co  be  fraely  assembled  I4,  25, 

(3)  Adding  newly  constructed  corrective  rational 
functions  to  the  basis  ^ These  rational 
functions  modify  the  smoothness  of  the  approxima¬ 
tions  at  vertices  (while  preserving  cl  continuity) 
but  permit  a  free  assembly  without  enforcing 
constraints.  A  method  has  been  devised  for  inte¬ 
grating  these  rational  functions  over  triangles 
directly  without  recourse  co  numerical  quadra- 
tive  16 . 


3.3.  Hierarchic  C  Family  of  Rectangular  Elements 


let  (3/3ni)  i  “  1,2,  be  Che  derivative  in  the  direction 
normal  co  side  i.  Then  in  order  for  the  approximation 
w  to  be  in  Cl  at  the  vertex’.che  following  constraint 
must  hold: 

32w  j2w 

—  cos  j>+  ~~  sin  p 

3s*  3slSnl 


Again  we  choose  as  nodal  variables  the  values  of 
the  approximation  u  at  vertices  and  higher  tangential 
derivatives  of  degree  j,  2  <_  j  <_  p  at  midside  nodes. 
First  observe  chat  che  polynomial 


Qj(C)  - 


j  2,  even 

j  ^  3,  odd 


3 " V  3 “y 

“  ~2  cos  *"  TsTfnT  3in  ♦  • 

3  s,  2  2 


Constraints  of  the  form  (3.1)  must  be  satisfied  at  all 
vertices  of  the  triangulacion. 

The  hierarchic  C^-  family  of  triangular  elemencs 
of  order  p  >  5,  uses  for  external  nodal  variables 
values  of  w,  its  first  derivatives  and  its  second 
tangential  and  mixed  (normal-tangencial)  derivatives 
1C  vertices,  and  derivatives  or  order  >  5  ic  midside 
nodes'-^  (See  Figure  10)  for  che  qulncic  hierarchic 
O  element).  For  each  6  <  )  ^  p,  the  jth  order 
cangenciai  derivatives  at  midsides,  and  a  mixed  jth 
order  derivative  (J-1  cangenciai  derivatives,  one  nor¬ 
mal  derivative)  at  midsides  are  used  to  enforce  C^- 
continuity. 


:u::<r:c  'I'-uiruut 
Hi'tnr 


hum  a/  :h,  dijoUcciwne  'wienon 
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3asls  functions  for  these  external  nodal  values  as  well 
as  for  internal  nodes  all  of  which  contain  a  factor 

4r®  iivan  ln  Several  procedures  are 

available  co  enforce  che  conscrainc  equations  (3.1): 


satisfies 

Qj(^l)  -  o  Q^(0)  -  0  i  -  2 . j-1 

Q^J)(0)  »  1  . 

S'ow,  consider  che  square  of  side  2  shown  in  Figure  11. 


Sasis  functions  corresponding  co  che  nodal 

variables  u(i)  1  -  1,2, 3, 4  are 

su(i)  ”  4  (L-°a-n>  \(3)  ■  xu+ou+n) 

MU(2)  ■  7  (1+«Hl-n)  su(4)  -  i(l-c)(i+n) 


and  it  is  easily  3een  that  chese  basis  functions  span 


iJ!3 
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che  sane  space  as  l,£,n,5n  i.e.  they  contain  Che  com¬ 
plete  linear 'polynomial.  Also  these  nodal  variables 
enforce  continuity  across  sides.  Mow  denoting  by 
(ij)  the  midpoint  of  side  ij,  basis  functions  corre¬ 
sponding  to  the  nodal  variables  urr(12),  u  (23)  , 
ur;(34),  unn(41)  are  ”  nn 

\  (12)  -  Q2CO(l-n)  S  (34)  -  Q2(-r)(l+n) 

55  55 

(3.3) 

S  (23)  -  (l+OQ,(ti)  tl  (41)  -  (1-5)CL(-t,) 

nn  4  nn  4 

and  chsse  basis  functions  added  to  the  ones  in  (3.2) 
span  the  same  space  as  l,5.n,52,5n,n2,52n,5n2,  i.e. 
chey  contain  a  complete  quadratic  and  they  enforce  Cu 
continuity  along  sides.  The  basis  functions  (3.2)  are 
taken  for  the  hierarchic  rectangular  cP  linear  element, 
and  those  in  (3.2)  and  (3.3)  for  the  quadratic  element. 

For  j  ^  3,  we  have 

Vj  (12)  *  Qj  (5)  (1_n)  Vj  (34)  ’  Qj  ("°  (1+n) 

^  (23)  -  (l+OQ^n)  Na.j(41>  *  (l-Ocy-n) 


as  the  basis  function  for  jth  order  tangential  deriva¬ 
tives  at  aidsides,  and  for  4  we  add  the  internal 
nodes 

(1  -  £2)  (1  -  q2)£j"4'V  t  -  0 . j— 4 . 


0 

Using  che  hierarchic  rectangular^C  elemencs,  it 
is  possible  to  construct  hierarchic  Cu  brick  elements, 
and  using  both  the  hierarchic  C°  rectangular  and.  CO 
triangular  elements  it  is  possible  to  construct  hier¬ 
archic  C°  triangular  prismatic  elements."8 
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5.  References 


These  basis  functions  span  Che  same  space  as  a  com¬ 
plete  polynomial  of  degree  and  the  two  monomials  of 
degree  j+1,  CJh,5h2.  Q 

Thus,  the  hierarchic  C  square  elemenc  of  degree 
P  JL  -  has  (1/2)  (p+i)  (p+2)+2  basis  functions,  two  more 
than  the  dimension  of  the  complete  polynomial  of 
degree  p.  The  cwo  extra  terms  correspond  to  5^n,5nP. 
3v  scaling  the  sides  of  che  3quare  the  elements  are 
easily  transformed  inco  rectangular  ones. 

3,i.  Hierarchic  C3  Solid  Elements 

Using  che  hierarchic  C®  triangular  and  rectangu¬ 
lar  elemencs,  we  can  construct  hierarchic  chree 


dimensional  elemencs  of  various  shapes.  For  exmaple, 
using  only  triangular  elements  we  construct  cetrahedral 
C®  elements  in  natural  (cetrahedral)  coordinates  (see 
for  example,  for  a  definition).  The  basis  func¬ 


tions  corresponding  to  vertices  of  all  triangular  faces 
of  the  tetrahedron  are  retained,  as  well  as  chose  cor¬ 


responding  to  sides  of  triangles  which  now  correspond 
to  edges  of  che  tetrahedron.  Those  corresponding  to 
incernal  nodes  of  triangles  now  correspond  to  face 


modes,  and  internal  modes  of  che  cecranedron  all  con¬ 


tain  the  factor  LjLpLjLq  i.e.  chey  vanish  on  all  faces 
of  the  tetrahedron  (3ee  Figure  12,  for  che  first  four 
hierarchic  cetrahedral  C®  elements) . 
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Comparative  Rates  of  h-  and  p- convergence  in  the 
Finite  Element  Analysis  of  a  Model  Bar 

_ Problem _ 

The  conventional  approach  to  finite  element  stress 
analysis  of  a  body  defined  by  a  polygonal  domain 
0  (in  two  dimensions)  is  to  triangulate  0  and  to 
seek  accuracy  by  letting  h,  the  maximum  diameter 
of  all  elements  in  the  triangulation,  tend  to 
zero.  This  approach,  called  h-convergence ,  has 
been  the  subject  of  intensive  investigation. 
Another  approach  which  is  being  developed  at 
Washington  University  is  to  ..fix  the  .triangulation 
of  and  to  let  p,  the  degree  of  the  complete, 
conf ormirfg ,  approximating  polynomial  over  each 
triangle,  tend  to  infinity.  Extensive  numerical 
tests  have  shown  that  the  second  approach,  called 
p- convergence ,  is  considerably  more’  accurate  than 
the  first,  even  in  problems  whose  solutions  have 
singularities  such  as  cracks  or  corners. 

In  order  to  illustrate  the  comparative  rates  of 
convergence,  a  model  (one-dimensional)  bar  prob¬ 
lem  is  studied.  Asymptotic  analysis  leads  to  ex¬ 
pressions  for  the  rates  of  convergence  in  the  two 
approaches,  when  the  solution  possesses  a  singu¬ 
larity  which  is  known  a  priori.  It  is  demon¬ 
strated  that  the  order  of  p-convergence  is  twice 
that  of  h-convergence,  provided  that  the  singular¬ 
ity  is  located  at  some  node  of  a  finite  element. 


